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PROGRAMME  DU  CONGRES 


Dimanche  21  mai 

18H30 

Cocktail  d'accueil  des  participants  a  I'Altea 
et  inscriptions 

Lundi  22  mai 

8H30 

Inscriptions  sur  le  lieu  du  Congres 
(Ecole  des  Mines  29  rue  Ponchardier) 

9H45 

Inauguration  du  congres  par  le  President  de 
rUniversite  et  le  Directeur  des  Mines 

lOH 

Information  sur  le  deroulement  du  Congres 

10H15-  11H05 
11H05-  11H20 
11H20-  11H50 
11H55-  12H25 

Conf.de  50'  GLIMM 

Pause  cafe 

Conf.de  30'  ANTONTSEV 

Conf.de  30'  CUSHMAN 

12H30  -  14H 

Repas 

15H30  -  17H 

Posters  1  a  13 

17H  -  17H50 

17H55  -  18H25 

19H 

Conf.de  50'  AVELLANEDA 

Conf.de  30'  GREENKORN 

Reception  a  la  Mairie  de  St  Etienne 

Mardi  23  mai 

9H  -  9H50 

9H55  -  10H25 

10H25  -  1 IH 

IIH  -11H50 

11H55-  12H25 

Conf.  de  50'  KNABNER 

Conf.  de  30'  ALLAIRE 

Pause  cafe  et  photo  sraphie  du  congres 
Conf.  de  50'  FASANO 

Conf.  de  30'  RUSSEL 

12H30  -  14H 

14H30 

15H30  -  17H 

Repas 

Conference  de  presse 

Posters  14  a  26 

17H  -  17H50 

17H55  -  18H25 

18H30  -  19H 

19H 

Conf.  de  50'  DOUGLAS 

Conf.  de  30'  LENORMAND 

Conf.  de  30'  GALLOUET 

Coktail  offert  par  I'Universite 

20H30 

Concert  offert  par  la  Mairie  choeur 
Contrechant  (salle  Aristide  Briand 
a  la  Mairie) 

Mercredi  24  mai 

9H  -  9H50 

9H55  -  10H25 

10H25-  IIH 

IIH  -  11H50 
11H55-  12H25 

Conf.  de  50'  EWING 

Conf.  de  30'  THOMAS 

Pause  cafe 

Conf.  de  50'  PANFILOV 

Conf.  de  30'  BADEA 

12H30  -  14H 

Repas 
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Mercredi  24  mai 

Apres-midi  libre : 

Prestations  proposees  aux  participants  (a  leur 
charge) : 

-  visite  du  Musee  d’Art  Modeme, 

-  visite  du  Musee  de  la  Mine, 

-  seance  au  planetarium 

-  golf, 

18h30 

Cocktail  et  presentation  de  I'economie  de  la 
region  stephanoise  a  la  Chambre  de 

Commerce  et  d'Industrie 

Jeudi  25  mai 

9H  -  9H50 

9H55  -  10H25 
10H25-  IIH 

IIH  -  11H50 
11H55-  12H25 

Conf.de  50'  WHEELER 

Conf.de  30’  RISEBRO 

Pause  cafe 

Conf.de  50'  MOLENAAR 

Conf.de  30'  YORTSOS 

12H30  -  14H 

Repas 

15H30-  17H 

Posters  2.'^  a-  4-^? 

17H  -  17H50 
17H55  -  18H25 
18H30  -  19H 

Conf.de  50'  VAZQUEZ 

Conf.de  30’  CHAVENT 

Conf.de  30’  GILBERT 

19H45 

Depart  des  cars  pour  le  Banquet  devant 

I'Altea  pour  le  Superflu  a  St  Romain  le  Puy 

Vendredi  26  mai 

9H  -  9H50 

9H55  -  10H25 
10H25-  IIH 

IIH  -  11H30 
11H35-  12H05 
12H10-  12H40 

Conf.  de  50'  ZHIKOV 

Conf.  de  30’  HORNUNG 

Pause  cafe 

Conf.  de  30’  AGANOVIC 

Conf.  de  30’  QUINT ARD 

Conf.  de  30’  GIPOULOUX 

12H45-  14H15 

Repas 
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PROGRAMME  SCIENTIFIQUE 
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50’  LECTURES 


AVELLANEDA  M.  Courant  Institute,  NEW  YORK  (U.S.A.) 

Exact  relations  between  electrical  measurements  and  permeability  :  a  porescale  analysis 

DOUGLAS  J.  PURDUE  University  (U.S.A.) 

Anomolous  diffusion  in  immiscible  displacement  resulting  from  fractal  permeabilities 


EWING  R.E  Texas  A  &  M  University,  COLLEGE  STATION  (U.S.A.) 
Mathematical  modelling  and  simulation  for  applications  of  fluid  flow  in  porous  media 


FASANO  A.  University  of  FIRENZE  (ITALY) 

Satured  flows  through  deformable  porous  media  with  mechanical  and  chemical  interactions 


GLIMM  J.  University  of  STONY  BROOK  (U.S.A.) 
Scaling  and  scale  up  for  flow  in  porous  media 


KNABNER  P.  University  of  ERLANGEN  (GERMANY) 

Finite  element  approximation  of  contaminant  transport  problems  in  porous  media 


MOLENAAR  J.  DELFT  University  of  Technology  (THE  NETHERLANDS) 
Density  driven  two-phase  flow  in  heterogeneous  porous  media 


PANFILOV  M.  University  OIL  &  GAZ  Research  Institute,  MOSCOU  (RUSSIA) 
Averaged  models  governing  tranfer  process  through  highly  heterogeneous  porous  media 


VAZQUEZ  L.  Universitad  Autonoma  de  MADRID  (SPAIN) 
Mathematical  model  for  two-phase  non  equilibrium  flows  in  porous  media 


WHEELER  M.  University  of  HOUSTON  (U.S.A.) 
non  parvenue 


ZHIKOV  V.V.  Pedagogical  University,  VLADIMIR  (RUSSIE) 
Homogenization  of  monotone  operations  in  perforated  domains 
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30’  LECTURES 


AGANOVIC  I.  University  of  ZAGREB  (CROATIE) 
Homogeneization  of  micropolar  flow  through  a  porous  medium 


ALLAIRE  G.  C.E.  A.  and  University  of  Paris  6,  PARIS  (FRANCE) 

Multi-scale  convergence  and  homogenization  in  porous  media  with  an  infmitnumber  of  length 
scales 


ANTONTSEV  S.N.  University  of  OVIEDO  (SPAIN) 

Filtration  of  immiscible  fluids  ;  effects  capillary  locking  and  hysteresis  self-propelling  under 
wave  impact 


BADEA  A.  Universite  Jean  Monnet,  ST  ETIENNE  (FRANCE) 
Homogenization  of  two-phase  flow  through  random  porous  media 


CHAVENT  G.  I.N.R.I.A.,  ROCQUENCOURT  (FRANCE) 

Methology  for  the  numerical  treatment  of  two  phase  flows  through  porous  media  with  two  rock- 
types 


CUSHMAN  J.H.  PURDUE  University  (U.S.A.) 

Comparaison  of  eulerian  nonlocal  models  to  lagrangian  models  of  chemical  transport  in 
heterogeneous  media 


ENE  H.I.  Institut  de  Mathematiques,  BUCAREST  (ROMANIA) 
On  the  homogenization  of  clays 


GALLOUET  T.  Ecole  Normale  Superieure,  LYON  (FRANCE) 

Schemas  de  volumes  finis  pour  les  problemes  d’ecoulement  en  milieu  poreux 


GILBERT  R.  University  of  Delaware,  NEWARK  (U.S'.A.) 

Boundary  layer  corrections  in  the  poroelastic  seabed  lying  under  a  shallow  ocean 


GIPOULOUX  O.  Universite  Jean  Monnet,  ST  ETIENNE  (FRANCE) 
Theorical  and  numerical  homogenization  of  polymer  flows 


GREENKORN  R.A.  PURDUE  University  (U.S.A.) 

An  examination  of  a  stochastic-nonlocal  theory  for  modeling  dispersion  in  scale-dependent 
porous  media 


HORNUNG  U.  University  B  W  M,  MUNCHEN  (GERMANY) 
A  model  of  multicomponent  diffusion  in  reactive  acid  soil 
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LENORMAND  R.  Institut  Fran9ais  du  Petrole,  RUEIL-MALMAISON  (FRANCE) 
Transport  equations  for  fluid  displacements  in  heterogeneous  porous  media  :  the  MHD  model 


QUINTARD  M.  E.N.S.A.M.,  BORDEAUX  (FRANCE) 

Calculation  of  porous  media  effective  properties  :  computational  problems  and  required  unit  cell 
feamres 


RUSSEL  T.F.  University  of  Colorado,  DENVER  (U.S.A.) 

Upscalling  of  dispervisity  in  modelling  of  solute  transport :  mathematical  theory  and  simulations 
of  laboratory  experiments 

THOMAS  J.M.  Universite  de  PAU  (FRANCE) 

Analysis  of  finite  volume  methods  and  application  to  reverse  simulation 


YORSTOS  Y.C.  University  of  Southern  California,  LOS  ANGELES  (U.S.A.) 
Studies  with  asymptotic  regimes  of  displacements  in  porous  media 


8 


POSTERS 


Monday  22th  of  May,  15  H  30 


1)  ON  A  MODEL  OF  FRACTURED  POROUS  MEDIA  p.45 

loana-Andreea  ENE  -  Jeannine  SAINT  JEAN  PAULIN 

2)  NON-LINEAR  DARCY'S  LAW  P -46 

Eduard  MARUSIC  PALOKA 

3)  INSTABILITY  THRESHOLDS  IN  MISCIBLE  FLUID  FLOWS  p.48 

D.  LOGGIA  -  N.  RAKOTOMALALA  -  D.  SALIN 

4)  HOMOGENEISATION  DES  EQUATIONS  DE  STOKES  ET  NAVIER-STOKES  DANS 
UN  MILIEU  POREUX  SANS  PROLONGEMENT  DE  LA  PRESSION  ET  ESTIMATION 
DES  CORRECTEURS 

EricBLAVIER  P-50 

5)  MIXED-HYBRID  FINITE  ELEMENT  APPROXIMATION  OF  THE  CHEMICAL 
SPACIES’  TRANSPORT  PROBLEM 

J.  MARYSKA  -  J.  MUZAK  P-52 

6)  DISPERSION  EN  MILIEUX  POREUX  PERIODIQUES  :  DE  LA  MODELISATION 

NUMERIQUE  A  L'EXPERIENCE  P •  54 

Sophie  DIDIERJEAN  -  Helio  AMARAL  SOUTO  -  Renaud  DELANNAY  -  Christian  MOYNE 

7)  CORRECTION  NON  LINEAIRE  DE  LA  LOI  DE  DARCY  :  ETUDE  NUMERIQUE  DE 
L'ECOULEMENT  BIDIMENSIONNEL  POUR  QUELQUES  ARRANGEMENTS 
PERIODIQUES 

Helio  AMARAL  SOUTO  -  Christian  MOYNE  P  •  5  6 

8)  MODELS  FOR  WATER  PERCOLATION  DURING  THE  PREPARATION  OF  ESPRESSO 
COFFE 

Gianni  BALDINI  P-58 

9)  MATHEMATICAL  MODELS  FOR  FREEZING  IN  POROUS  MEDIA 

Federico  TALAMUCCI  P-59 

10)  DIFFUSION  AND  DISSOLUTION  IN  REACTIVE  POROUS  MEDIUM  ; 
MATHEMATICAL  MODELING  AND  NUMERICAL  SIMULATIONS 

E.  MAISSE  -  P.  MOSZKOWICZ  -  J.  POUSIN  -  F.  SANCHEZ  p.61 

1 1)  LATTICE  GAS  SIMULATION  OF  VISCOUS  FINGERING  FOR  MISCIBLE  FLUIDS 

N.  RAKOTOMALALA  -  D.  SALIN  -  P.  WATZKY  p .  64 

12)  NONLINEAR  DIFFUSION  IN  AN  INHOMOGENEOUS  AQUIFER 

M.  GUEDDA  -  D.  HILHORST  -  M.  PELETIER  P-66 

13)  SPATIALISATIONS  QUANTITATIVE  ET  QUALITATIVE  EN  MILIEUX  POREUX 

APPLICATION  A  L'ESTIMATION  DES  RISQUES  DANS  LA  DEPOLLUTION  DES  SOLS 
W.  ANKER  -  D.  GRAILLOT  -  J.  BOURGOIS  -  M.  ZELFANI  p.68 
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Tuesday  23th  of  May,  15  H  30 


14)  WATER  RETENTION  CHARACTERISTICS  FOR  FRACTAL  SOIL 

Nigel  BIRD  p.70 

15)  CONTRIBUTION  OF  IMAGE  ANALYSIS  TO  ACCURATE  NUMERICAL  FLOW 
MODELLING 

Y.  ANGUY  -  D.  BERNARD  -  R.  EHRLICH  p.7 1 

16)  A  THREE  DIMENSIONAL  MODEL  FOR  FLOW  AND  TRANSPORT  IN  SATURATED 
POROUS  MEDIA  APPLICATION  TO  SALT  INTRUSION  INTO  COASTAL  AQUIFER 

F.  JACOB  -  J.M.  CROLET  -  P.  LESAINT  -  J.  MANIA  p.73 

17)  USE  OF  THE  FOURIER-LAPLACE  TRANSFORMATION  AND  OF 
DIAGRAMMATICAL  METHODS  TO  INTERPRET  PUMPING  TEST  IN 
HETEROGENEOUS  RESERVOIRS 

Benoit  NOETINGER  -  Y.  GAUTIER  p.75 

18)  COMPRESSIBLE  FLOW  IN  POROUS  MEDIA.  CHECKING  OF  A  MACROSCOPIC 
MODEL  BY  MEANS  OF  PORE-LEVEL  SIMULATIONS 

A.  BOUHOUCH  -  M.  PRAT  -  H.  BOISSON  p.76 

19)  MATHEMATICAL  MODELLING  OF  DIAGENETIC  PROCESSES  IN  SEDIMENTARY 
BASINS 

Astrid  HOLSTAD  p .  77 

20)  MEDIAS  AXIS  ANALYSIS  OF  THREE  DIMENSIONAL  TOMOGRAPHIC  IMAGES  OF 
DRILL  CORE  SAMPLES 

W.  BENT  LINDQUIST  -  Sang-Moon  LEE  p.79 

21)  ETUDE  NUMERIQUE  DE  LA  METHODE  DE  RECUPERATION  SECOND  AIRE  DU 
PETROLE 

MazenSAAD  p.81 

22)  SCALING-UP  PERMEABILITY  IN  THE  NEAR  WELL  REGIONS 

Y.  DING  p.83 

23)  ON  THE  DAM  PROBLEM  WITH  LEAKY  BOUNDARY  CONDITIONS  AND  LE4ER 
DARCY'S  LAW 

A.  LYAGHFOURI  p.85 

24)  NUMERICAL  MODELLING  OF  CLAY  DRYING 

E.F.  KAASSCHIETER  p.86 

25)  THE  PLANE  POTENTIAL  FLOW  THROUGH  THE  INHOMOGENEOUS  POROUS 
LINE 

Ivan  KEGLEVIC  p.87 

26)  IMPLEMENTING  A  NON-LOCAL  FOR  FLOW  AND  TRANSPORT  THROUGH 
POROUS  MEDIA 

John  F.  PETERS  -  Stacy  E.  HOWINGTON  p.89 
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Thursday  25th  of  May,  15  H  30 

27)  EFFECTIVE  PERMEABILITY  OF  STRONGLY  HETEROGENEOUS  POROUS  MEDIA 

A.DeWIT  P-91 

28)  DIFFERENT  PRESSURE  GRIDS  IN  HETEROGENEOUS  POROUS  MEDIA 

Sophie  VERDIERE  -  Dominique  GUERILLOT  p . 93 

29)  CONTRIBUTION  A  L'ETUDE  DU  COMPORTEMENT  DE  LA  SOLUTION 
DEQUATIONS  NON  LINEAIRES  DE  DIFFUSION-CONVECTION 

J.B.  BETBEDER  -  L.  LEVI  -  A.  PLOUVIER  -  G.  VALLET  p.95 

30)  SUR  LA  CORRECTION  NON-LFNEAIRE  DE  LA  LOI  DE  DARCY  POUR  LES 

ECOULEMENTS  EN  MILIEUX  POREUX  A  FAIBLE  NOMBRE  DE  REYNOLDS 
Mouaouia  FIRDAOUSS  -  Jean  Luc  GUERMOND  p.97 

31)  THE  CAUCHY  PROBLEM  FOR  HYPERBOLIC  CONVERSATION  LAWS  WITH 
THREE  EQUATIONS 

Yun-guang  LU  -  Christian  KLINGENBERG  p.  1 0 1 


32)  MODELLING  OF  A  HORIZONTAL  AND  A  VERTICAL  FRACTURED  WELL 

B.  CVETKOVIC  -  G.  HALVORSEN  -  E.  LOW  p.  102 

33)  RESOLUTION  NUMERIQUE  D'UN  PROBLEME  DE  RESTAURATION  BIOLOGIQUE 
EN  MILIEUX  POREUX 

Ch.  H.  BRUNEAU  -  P.  FABRIE  -  P.  RASETARINERA  p.  103 

34) ASYMPTOTIC  BEHAVIOUR  OF  SOLUTIONS  OF  NEUMANN  PROBLEMS  FOR 
NONLINEAR  NONVARIATION  ELLIPTIC  EQUATIONS  IN  DOMAINS  WITH 
COMPOSITE  STRUCTURE 

Alexander  A.  KOVALEVSKY  p.  1 05 

35)  STOKAGE  DE  GAZ  REACTIF  DANS  UN  MILIEU  POREUX  MODELISATION 
MATHEMATIQUE  ET  NUMERIQUE 

Victor  DUVAL  -  Hamid  GHIDOUCHE  -  Claude  B ASDEVANT  p .  1 06 

36)  MACROSCOPIC  PERMEABILITY  OE  THE  SYSTEM  OF  THIN  FISSURES  FILLED 

BY  MATERIAL  WITH  RANDOM  PERMEABILIIY  TENSOR  p .  1 07 

G.P.  PANASENKO 

37)  APPLICATION  OF  CONVECTIVE  TRANSFERS  ON  GEOLOGY 

B.  GERARD  -  JJ.  ROYER  -  C.  LE  C.ARLIER  DE  VESLUD  p.  109 

38)  THE  EXPLICIT  SOLUTION  OF  A  FREE  BOUNDARY  PROBLEM  FOR  A 

NONLINEAR  ABSORPTION  MODEL  OF  MIXED  SATURATED-UNSATURATED  FLOW 
A.C.  BRIOZZO  -  D.A.  TARZIA  P- 1 1 1 

39)  ON  SOME  NEW  MODELS  OF  NONHOMOGENEOUS  FLUID  FILTRATION  IN 
POROUS  MEDIA 

V.N.  MONAKHOV  P-112 

40)  MULTISCALE  MODELS  :  A  TOOL  TO  DESCRIBE  THE  POROSIIY  OF  CEMENT- 
BASED  MATERIALS  AND  TO  PREDICT  THEIR  TRANSPORT  PROPERTIES 

J.F.  DAIAN  -  D.  QUENARD  P-113 
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CONFERENCES 
DE  50 
MINUTES 


EXACT  RELATIONS  BETWEEN  ELECTRICAL  MEASUREMENTS 
AND  PERMEABILITY:  A  PORE-SCALE  ANALYSIS 

Marco  Avellaneda 

Courant  Institute,  New  York  University 

Electrical  measurements  have  been  used  for  many  years  by  petroleum  engineers 
to  characterize  the  formation  and  fluid  conductivity  of  porous  media.  In  this  talk, 
relations  between  electrical  measurements  and  permeability  are  discussed  and  an¬ 
alyzed  through  simple  models  and  computer  simulations  of  Stokes  flow  through 
capillaries.  It  is  shown  that  certain  estimators  based  on  electric  measurements  pro¬ 
posed  by  Johnson  et  al.  (Jr.  Fluid  Mech.,  176,  379  (1986))  and  Zhou  and  Sheng 
(Phys.  Rev.  Lett.  57,  2565.  (1986))  are,  generically,  good  predictors  of  the  effective 
pore-radius  and  the  Darcy  constant.  However,  the  accuracy  of  such  estimators  may 
deteriorate  if  the  porous  medium  has  a  wide  range  of  (active)  pore  sizes  or  due  to 
the  asperity  of  pore  walls. 


Anomolous  Diffusion  in  Immiscible 
Displacement  Resulting  from  Fractal 
Permeabilities 

Jim  Douglas,  Jr.  *  Felipe  Pereira*  Frederico  Furtado^ 

March  15,  1995 


Abstract 

The  analysis,  through  high  resolution  numerical  simulations,  of 
mixing  lengths  in  water-oil  fronts  in  heterogeneous  porous  media  can 
be  accomplished  through  efficient  computations  on  fast,  multi-processor 
computers.  This  paper  is  devoted  to  the  description  of  an  algorithm 
for  such  a  computation  and  to  the  presentation  of  numerical  results 
that  lead  to  the  characterization  of  anomalous  diffusion  arising  from 
fractal  permeabilities.  It  is  shown  that  the  exponent  7  describing  the 
scaling  behavior  (for  large  time)  of  the  mixing  region  lies  in  the  range 
5  <  7  <  1-  For  physically  typical  heterogeneity  strengths  and  mo- 
bihty  ratios  occurring  in  waterflooding  processes,  7  =  1;  thus,  typical 
waterflooding  displays  non-Fickean  behavior,  while  the  behavior  tends 
to  Fickean  as  the  mobility  ratio  tends  to  one. 

The  computational  algorithm  is  based  on  a  domain-decomposition 
technique  apphed  to  a  mixed  finite  element  approximation  of  the  non- 
bnear  equations  for  two-phase,  immiscible,  incompressible  flow  ex¬ 
pressed  in  an  Eulerian  form  and  is  implemented  to  be  portable  to 
several  state-of-the-art  parallel  architectures. 


’Center  for  Applied  Mathematics,  Purdue  University,  West  Lafayette,  IN,  USA 
^Department  of  Mathematics,  UNICAMP,  Campinas,  Brazil 


MATHEMATICAL  MODELING  AND  SIMULATION 
FOR  APPLICATIONS  OF  FLUID  FLOW 
IN  POROUS  MEDIA 


Richard  E.  Ewing 

Institute  for  Scientific  Computation 
Texas  A&M  University 


Abstract 


Understanding  the  fate  and  transport  of  contaminants  to  determine  wa¬ 
ter  quality  and  to  develop  remediation  strategies  or  optimizing  the  recovery  of 
hydrocarbons  in  petroleum  applications  each  require  the  ability  to  model  multi¬ 
phase  flow  in  heterogeneous  three-dimensional  reservoirs.  Model  equations  and 
corresponding  parameters  must  be  determined  at  the  appropriate  length  scales 
to  describe  the  scaled  physics  of  flow.  Effective  simulators  require  accurate  nu¬ 
merical  methods  on  general  geometries.  Use  of  mixed  finite  element  methods 
and  local  grid  refinement  w’ill  be  discussed.  Example  calculations  for  field  sim¬ 
ulations  in  aquifers  or  reservoirs  with  complex  boundaries  will  be  presented. 
Parallelization  of  the  codes  will  also  be  discussed. 
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Saturated  flows  through  deformable  porous  media 
with  mechanical  and  chemical  interactions 

Antonio  FASANO 
University  of  Firenze,  ITALY 

We  will  examine  a  variety  of  filtration  processes  in  which  there  is  a  mutual 
interaction  between  the  flow  and  the  porous  medium.  Typical  cases  are  the  dissolution 
of  substances  from  the  porous  matrix  to  the  flow,  the  removal  of  small  particles  from  the 
skeleton  and  their  transport  by  the  flow,  the  deformation  of  the  medium  induced  by  the 
flow  with  the  consequent  change  of  the  porosity  and  of  the  hydraulic  conductivity.  In 
the  above  framework  many  different  situations  can  arise.  For  instance  the  deformations 
can  be  reversible  or  not,  they  can  be  instantaneous  or  can  obey  some  relaxation  kinetics. 

Also  we  can  have  a  wetting  front  (i.e.  a  free  boundary  problem)  or  other  types  of  free 
boundaries  due  to  the  generation  of  internal  inhomogeneities  because  of  the  transport 
of  solid  components.  The  variety  of  the  physical  situations  is  obviously  reflected  in  the 
diversity  of  the  mathematical  structure. 


Scaling  and  Scale  up  for  Flow  in  Porous  Media 

James  GLIMM 

University  at  Stony  Brook,  Stony  Brook  NY  11794-3600 

Geological  variation  occurs  across  all  length  scales.  Both  fine  and  coarse  scale 
variability  produce  important  effects  upon  the  coarse  scale  fluid  flow,  which  is  of  direct 
engineering  interest.  The  (scale  up)  problem  of  predicting  the  influence  of  fine  scale 
features  on  coarse  scale  flow  is  further  complicated  by  the  fact  that  the  fine  scale 
variation  is  only  specified  stochastically,  as  a  random  field. 

As  a  simple  model  of  multi-length  scale  geological  variation,  we  use  fractal  random 
fields.  For  the  case  of  linear  transport  (flow  of  a  passive  scalar  concentration),  results 
of  the  author  and  colleagues  will  be  presented,  based  on  direct  simulation  and  ordinary 
and  renormalized  perturbation  theory.  Comparison  to  results  of  others  will  be  included. 
Results  of  the  author  and  others  concerning  nonlinear  stochastic  transport  (which  is  a 
considerably  more  difficult  problem)  will  also  be  presented. 


Finite  element  approximation  of  contaminant 
transport  problems  in  porous  media 


Peter  Knabner 

University  of  Erlangen-Niirnberg 
Institute  for  Applied  Mathematics 
Martensstr.  3,  91058  Erlangen 
Germany 

Models  for  the  transport  of  reactive  solutes  in  porous  media  as  the  adsorption- 
diffusion-advection  model  or  more  involved  multi-component  models  exhibit 
various  aspects  which  distinguish  them  from  smooth  parabolic  equations  or 
systems; 

-  non-smooth  nonlinearities,  as  e.q.  the  Freundlich  isotherm, 

-  non-equilibrium, 

-  advection  dominance. 

I.e.,  emerging  mathematical  problems  are  nonlinear  parabolic  systems  with 
various  degenerations: 

-  parabolic-hyperbolic  locally 

(as  in  the  porous  medium  equation), 

-  parabolic-ordinary  differential  equation, 

-  parabolic  close  to  hyperbolic  globally 
(as  singular  limit). 

In  the  first  part  of  the  lecture  we  will  discuss  order  of  convergence  results 
for  the  conformal  finite  element  method  taking  the  first  two  phenomena  into 
account.  We  allow  for  numerical  quadrature  and  time  discretization  and 


elucidate  the  role  of  regularization  and  kinetic  relaxation  of  equilibrium  re¬ 
actions.  If  in  addition  the  third  aspect  is  prominent,  the  conformal  finite  ele¬ 
ment  approach  is  no  longer  appropriate  and  we  select  the  Lagrange-Gelerkin 
method  (or  modified  method  of  characteristics).  The  main  problem  con¬ 
sists  in  a  self-consistent  definition  of  the  chracteristics  to  insure  an  accurate 
tracking  of  fronts  and  shocks  also  for  large  time  steps.  The  evolving  “fully 
implicit  scheme”  is  discussed  and  preliminary  order  of  convergence  estimates 
indicated. 

-  joint  work  with  John  Barrett  (London)  and  Holger  Kappmeier  (Erlan¬ 
gen). 
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Density  driven  two-phase  flow  in  heterogeneous  porous  media 


J.  Molenaar 

TWI,  Delft  University  of  Technology, 
p.o.box  5031.  2600  GA  Delft,  The  Netherlands 


The  contamination  of  groundwater  by  Dense  Non-Aqueous  Phase  Liquids  is  increasingly  stud¬ 
ied  by  hydrologists.  A  DNAPL  exists  as  a  separate  phase  in  the  subsurface,  and  does  not 
mix  with  the  water  on  time  scales  of  practical  interest.  Therefore  the  spreading  of  a  DNAPL 
in  the  soil  can  be  described  by  an  immiscible,  incompressible  two-phcise  flow  model.  In  the 
absence  of  capillary  forces  the  ID  flow  through  a  homogeneous  porous  medium  is  modeled  by 
the  classical  Buckley-Leverett  equation. 

In  this  presentation  we  consider  the  flow  of  a  DNAPL  in  a  heterogeneous  porous  medium 
taking  capillarity  into  account.  It  is  well  known  that  at  heterogeneity  a  certain  capillary  pres¬ 
sure  must  built  up,  before  the  DNAPL  can  enter  a  low  permeability  region.  The  minimum 
pressure  needed,  is  called  the  displacement  or  threshold  pressure.  A  nonzero  displacement 
pressure  in  the  capillary  pressure  model  (e.g.  the  Brooks-Corey  model)  has  interesting  math¬ 
ematical  consequences.  By  a  regularization  technique  we  derive  conditions  to  match  the 
solution  on  both  sides  of  an  interface,  where  the  soil  properties  are  discontinuous.  There  are 
two  conditions:  a  flux  condition  and  an  extended  capillary  pressure  condition.  To  illustrate 
the  implications  of  these  interface  conditions  we  consider  the  stationary  ID  problem.  With¬ 
out  gravity  this  problem  has  a  unique  steady-state  solution,  however  with  gravity  there  are 
multiple  steady-state  solutions.  Moreover  neither  the  DNAPL  saturation,  nor  the  capillary 
pressure  need  to  be  continuous  at  the  heterogeneity.  This  means  that  the  usual  assumptions 
for  standard  numerical  simulators  are  in  general  not  valid. 

We  present  a  numerical  algorithm  that  is  developed  specially  to  deal  with  this  situation: 
at  a  heterogeneity  we  assume  continuity  of  flux  and  the  extended  interface  condition  for 
the  capillary  pressure.  To  validate  the  numerical  method  we  consider  two  problems  with 
heterogeneities  that  can  be  treated  analytically:  the  ID  steady  state  situation  with  gravity, 
and  a  2D  time-dependent  problem  without  gravity.  Finally  we  will  show  some  computational 
results  for  cases  which  are  of  more  practical  relevance. 
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Averaged  Models  Governing  Transfer  Processes 
Through  Highly  Heterogeneous  Porous  Media 

M.  PANFILOV 
Moscow,  RUSSIA 

Macroscopic  models  are  proposed  to  describe  convection-diffusion  transfer  through 
multiply  heterogeneous  media.  The  heterogeneity  is  given  by  three  periodic  fields  of 
porosity,  conductivity  and  dispersion  parameter.  The  high  heterogeneity  corresponds 
to  the  case,  when  the  oscillation  amplitude  of  one  of  fields  is  large.  Dual  porosity 
medium  presents  one  example  of  examined  systems.  The  convection  velocity  is  derived 
from  the  pressure  distribution  described  by  parabolic  equation. 

All  media  are  shown  to  be  nonuniform.  They  include  15  classes,  with  different 
averaged  models  and  internal  flow  structures.  Media  when  the  convection  predominates 
in  blocks  present  essential  interest.  They  can  be  decomposed  into  three  sub-classes 
depending  on  convection  flow  type.  Closed  results  were  obtained  for  the  through-type 
flow  in  blocks  or  source-type  flow  using  homogenization  method.  Averaged  models 
contain  the  exchange  term  that  is  no  symmetric  in  bloc  and  in  fracture.  Homogenization 
results  depend  upon  the  direction  of  pressure  variation  in  time.  In  case  when  the  time 
derivative  of  pressure  changes  its  sign,  the  averaged  model  consists  of  equations  with 
delaying  argument.  Example  of  its  solution  is  constructed.  It  corresponds  to  the  special 
inverse  method  of  oil  well  investigation  by  means  of  tracer  injection. 

Media  with  through-source  flow  type  in  blocs  was  examined  by  the  way  of 
simulations  on  the  network  model  of  porous  media.  The  analysis  of  exchange  process  is 
shown.  Some  results  for  nonlinear  case  were  presented  also.  This  system  corresponds  to 
both  two-phase  flow  or  one-phase  concentration  transfer,  when  fluid  viscosity  depends 
on  concentration.  The  case  of  two-phase  flow  was  examined  without  a  contribution  of 
high  heterogeneity.  Generalization  of  Darcy’s  law  including  the  capillary  dynamics  is 
derived.  It  is  shown  by  averaging  of  numerical  simulations,  that  new  type  of  capillary 
pressure  function  has  a  form  of  non  monotone  curves. 


Title:  Mathematical  model  for  two-phase 
non-equilibrium  flows  in  porous  media 

Juan  Luis  Vazquez,  Dpto.  de  Matematicas, 

Univ.  Autonoma  de  Madrid,  28049  Madrid,  Spain 

Abstract 

We  study  a  model  for  the  flow  of  two  immiscible  fluids  in  a  homogeneous  and  isotropic 
porous  medium,  based  on  the  Muskat-Leverett  model,  classical  in  the  oil  science,  with 
the  non-equilibrium  effects  introduced  by  Barenblatt  in  describing  the  process  of  water-oil 
displacement. 

It  leads  to  the  following  system  for  the  saturation  s  G  (0, 1)  and  the  effective  saturation 
a: 


eA(f){a)  =  a  —  s, 

TSt  +  S  =  (7, 

where  e  and  r  are  positive  constants  and  #  is  a  monotone  constitutive  function. 

We  pose  the  problem  for  x  >  0  and  t  >  0  with  intial  data  s(x,  0)  >  0  and  boundary  data 
(7(0,  f)  =  1.  Our  main  results  concern  the  existence  and  properties  of  the  free  boundary 
or  leading  front  of  the  set  {s(.T.f)  >  0}  and  the  large-time  behaviour  of  the  solution  and 
free  boundary. 

This  is  joint  work  with  G.I.  Barenblatt  (Urbana,  USA),  and  .Jestis  Garcia  Azorero  and 
Arturo  de  Pablo,  from  Madrid.  Previous  analysis  is  due  to  G.I.  Barenblatt  and  A.  Gilman. 


Homogenization  of  monotone  operators  in 
perforated  domains 

V.V.Zhikov 
January  10,1995 


Vladimir  Pedagogical  University 
Department  Mathematics 
600024,  Vladimir,  Russia 


Let  Q  be  a  connected  periodic  domain  in  (0.1)  ^ 

cell  of  periodicity.  We  consider  the  problem 

(  -div(a(£-ii,  Vup  +  b{e-^x,  u^))  =  /  in  fi  n  Q., 
i  a(£-ir,  Vm)  •  nlnnag,  =  0  (Neumann’s  condition), 
i  u'lannQ.  =  0  (Dirichlet’s  condition). 

Here  D  is  a  bounded  domain  in  IR^  and  a(i,, 0.  verifies  the  following 

structure  conditions 

1)  a{-,0<b{-,s)  are  D-periodic  and  Lebesgue  measurable  on  Q; 

2)  for  a.e.  x  eQ  and  every  ^i,6  £ 

(a(x,4i)  -  a(x,^2))  •  (6  -  ^2)  >  «  >0.  \  in  the  case  p>‘2 

la(r,4l)  -  a(x,^2)l  <  ““Hi  +  l^il^  +  l^^l)  )16-6lJ 


(a(x,6)  -  a(2;.6))  •  (^1  -  ^2)  >  “I'ei  -  ^ 

|a(x,fi)-a(x,6)l  <  “  -^2r 


in  the  case  1  <  p  <  2 


3)  a(x,0)  =  6(^.0)  =  0 

4)  for  a.e.  x  eQ  and  for  every  si,  S2  £  IR^ 

(6(x,  si)  -  6(x,  S2))  •  (si  —  •S2)  >  gHi  -  •52!^ .  /?  >  0.  \  in  the  case  g  >  2 

16(x,si)  -  K2:.«2)1  < /?“Hl  + l^il^  )lsi-S2li 

(6(x, si)  —  6(x, S2))  ■  (si  —  S2)  >  /ll^i  —  S2I  (l^ij  +  1^2!)  1  (.3.30  1  <  g  <  2 

16(x,Si)-6(x,S2)|<rHsi-S2l'  J 


Note  that  there  is  not  connection  between  the  exponents  p  and  q. 

The  homogenized  operator  gq  :  is  defined  for  every  ^  G  IR^ 

by 

^o(^)  =  [  a{y,(  +  Dw{y))dy, 
anQ 

where  w  is  the  unique  solution  to  the  problem 

f  I  “(?/)  +  D“-’)  •  ^<t>dy  =  0  for  every  4>  £ 

[wew^fAQ)- 


Let  also  set 


1  if  a;  G  <3 
0  if  X  G  IR^'^  \  <3 


bo{s)  =  y  x(x)6(x,5)dx,  0  =  \nr\Q\ 

a 

Theorem.  If  f  £  q'  —  then  we  have 

lim  /  =  0, 

5-0  y  ' 

OnQ, 

y(£~^x)a(c“^x,  Vu“)  ao(Vu°)  weakly  inL^{Il) 
where  is  the  unique  solution  to  homogenized  equation 

(  —div(ao{'Vu^)  +  bo{u°))  —  9f  in 
I  u°|aQ  =  0. 

The  case  of  <3  =  IR'^  see  Tartar  [1]  and  Fusco,  Moscariello  [2].  The  case 
of  variational  problem  in  perforated  domains  see  Zhikov  [3]. 
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Homogenization  of  Micropolar  Flow  through  a  Porous  Medium 

I.  Aganovic,  Zagreb 

Homogenization  of  an  incompressible  micropolar  fluid  flow  through  a  pe¬ 
riodic  porous  medium  is  discussed.  The  homogenized  problem  consists  of  the 
corrected  Darcy  law  and  divergence  free  condition  for  macroscopic  velocity. 
Under  appropriate  scaling  the  corresponding  convergence  result  is  proved. 

The  presented  results  are  a  part  of  the  joint  research  of  I.  Aganovic  and 
Z.  Tutek. 


MULTI-SCALE  CONVERGENCE  AND  HOMOGENIZATION 
IN  POROUS  MEDIA 

WITH  AN  INFINITE  NUMBER  OF  LENGTH  SCALES 

Gregoire  ALLAIRE 
Commissariat  a  I’Energie  Atomique 
DRN/DMT/SERMA,  C.E.  Saclay,  F-91191  GIF  sur  YVETTE 

k, 

Laboratoire  d’Analyse  Numerique,  Universite  Paris  6 


In  collaboration  with  M.  Briane,  we  generalized  the  notion  of  two-scale  con¬ 
vergence  to  the  case  of  multiple  separated  scales  of  periodic  oscillations.  We 
introduced  a  multi-scale  convergence  method  for  the  so-called  reiterated  ho¬ 
mogenization  of  partial  differential  equations  with  oscillating  coefficients.  This 
new  method  is  applied  to  a  diffusion  problem  in  a  porous  media  modeled  by  an 
infinite  number  of  periodic  scales  of  heterogeneities. 


Filtration  of  Immiscible  Fluids:  Effects  of 
Capillary  Locking  and  Hysteresis  Self-Propelling 

under  Wave  Impact 


Antontsev  S.N. 

Lavrentyev  Institute  of  Hydrodynamics,  Novosibirsk,  Russia 
FICYT,  University  of  Oviedo,  Spain 

Penkovskii  V.I. 

Lavrentyev  Institute  of  Hydrodynamics,  Novosibirsk,  Russia 


There  are  studied  qualitative  properties  of  solutions  of  the  classical  model  of  filtration 
of  two  immiscible  incompressible  liquids 

Vi  =  -Ko{x)ki{si)/ fj.iV{pi  +  pigh),  i  =  1,  2, 

+  div(p,rT)  =  0,  (z  =  1, 2,  5i  +  52  =  1) 

P2-P\^  Pc{si,x). 

It  is  shown  that  the  domain  of  filtration  flow  may  contain  ’’dead"  cores  of  two  types  where 
either  U  =  0  or  172  =  0.  The  former  happens  in  the  presence  of  the  limiting  values  of  the 
saturation  5°,  which  make  the  relative  phase  permabilities  zero,  ki(s°)  =  0.  Besides, 
the  saturation  distribution  in  the  dead  core  is  constant  (either  s  =  s°  or  s  =  5°). 

The  second  type  is  caused  by  the  action  of  the  capillary  locking  effect.  In  the  cores  of 
this  type  the  saturation  needn’t  be  constant  but  so  is  the  presssure  in  the  corresponding 
phase,  (pi  =  const.).  It  is  shown  that  among  with  the  well-known  effect  of  the  ’’boundary 
locking”  of  the  wetting  phase  an  analogous  phenomenon  can  occure  inside  the  flow  domain. 

There  are  considered  peculiarities  of  the  filtration  flow  of  two  immiscible  incompress¬ 
ible  fluids  in  mixed-wetted  porous  media  where  the  capillary  pressure  function  Pc{si)  is 
of  the  hysteresis  character.  Besides,  the  effect  of  the  wave  change  of  the  pressure  at  the 
production  and  injection  wells  is  discussed. 

There  are  presented  some  results  of  numerical  and  physical  experiments, 
where  u,-,  p,,  p,-  and  s,-  are  the  velocity  ,  tensity,  pressure,  saturation  of  the  phase  i, 
i  =  1,2,  m{x)  is  porosity  of  the  medium,  Ko{x)  is  the  tensor  of  absolute  permeability, 
ki{s)  —  koi/pi  ,  Pi  is  the  viscosity,  ^:oi(5i)  >  0  (^^{(O)  —  0,A:o,(5,)  increasing  in  s,)  is  the 
relative  permeability,  g  is  the  acceleration  due  to  gravity,  h  is  the  distance  to  a  fixed 

horizontal  reference  plane,  qi{x,t),  {qi  +  q2  =  0),  is  the  injected  mass  of  the  phase  f,  Pc  is 
the  capillary  pressure  (increasing  in  5i  [dp^lds  <  0)  and  p{x,Si)  =  0). 
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HOMOGENIZATION  OF  TWO-PHASE  FLOW  THROUGH 
RANDOM  POROUS  MEDIA 


by  Anca  BADEA^  and  Alain  BOURGEATl 


We  are  interested  to  apply  the  theory  of  homogenization  for  random  coefficients, 
developped  mainly  for  linear  operators,  to  the  modelisation  of  incompressible  two-phase 
flow  in  randomly  heterogeneous  porous  media. 

After  giving  the  equations  of  two-phase  flow  in  porous  media,  we  define  the  problem 
of  homogeneization  (or  scaling  up)  we  are  considering. 

In  a  second  part,  we  recall  the  theoretical  results  on  this  problem  as  obtained  in  [1] 
and  in  [2]. 

Due  to  the  difficulty  to  use  these  theoretical  results  as  a  mean  of  computing  the 
effective  coefficients  (permeability  and  porosity),  in  the  ergodic  case,  we  design  a  way  of 
doing  this  scaling-up  numerically  by  using  the  notion  of  Volume  Averaging.  We  prove  that 
this  method  gives  back  the  theoretical  effective  coefficients. 

In  the  last  part  we  present  numerical  tests.  The  first  test  is  for  a  one-phase  flow 
problem  and  shows  the  convergence  of  the  method.  The  second  one  is  for  two-phase  flow 
and  shows  interaction  between  fingering,  heterogeneities  and  scaling-up. 
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Mixed  finite  elements  and  cell-centered  finite  volumes 
for  two-phase  flow  in  porous  media  with  two  rock-types 

Guy  CHAVENT 

INRIA,  Rocquencourt,  FRANCE 

Finite  Volumes  methods  are  very  popular, as  they  are  very  close  to  the  physics 
and  allow  the  implementation  of  complicated  physical  laws.  However  their  design  may 
become  delicate  in  complex  geometries.  On  the  other  hand,  Finite  Elements  methods 
are  less  close  to  the  physics,  but  they  come  with  a  good  mathematical  background,  and 
their  analysis  is  quite  advanced  for  structured  as  well  unstructured  meshes. 

We  show  how  to  combine  the  advantages  of  finite  volumes  and  finite  elements 
methods  into  a  numerical  procedure  for  the  resolution  of  two-phase  incompressible  flows 
through  a  porous  medium  with  two  different  rock- types,  by  a  joint  use  of  a  mixed-hybrid 
finite  elements  method  and  of  a  Godunov  method. 
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Comparison  of  Eulerian  Nonlocal  Models  to  Lagrangian  Models  of 
Chemical  Transport  in  Heterogeneous  Media 

John  H.  Cushman.  BiU  X.  Hu,  and  Fei-Wen  Deng 

Abstract 

When  developing  transport  theory  for  heterogeneous  porous  media  there  are 
basically  two  frame  works  one  may  adopt;  Eulerian  and  Lagrangian.  The  standard  starting 
point  for  both  types  of  theories  is  a  “Darcy-scale”  CDE.  In  the  particle  tracking 
(Lagrangian)  framework,  the  Darcy-scale  dispersive  term  in  the  CDE  is  neglected  so  that 
the  particles  move  along  streamlines.  In  this  approach  one  does  not  obtain  the  mean 
concentration,  but  rather  various  order  mean  spatial  moments.  Using  the  Eulerian 
perspective,  there  is  no  need  to  neglect  Darcy-scale  dispersion  as  particles  are  not  tracked 
along  streamlines,  but  rather  they  are  viewed  from  a  fixed  position.  Thus  as  currently 
employed  the  Eulerian  framework  which  explicitly  obtains  the  mean  concentration  is  more 
general  than  the  Lagrangian  models.  We  show  that  if  the  local-scale  dispersion  is 
neglected  in  the  Eulerian  models  and  if  these  models  are  not  localized  for  reactive 
transport,  then  the  mean  spatial  moments  agree  with  those  obtained  through  the 
Lagrangian  picture.  However,  if  the  local-scale  dispersive  flux  is  kept  in  the  Eulerian 
model,  then  the  Eulerian  moments  can  be  substantially  different  than  the  Lagrangian 
moments. 
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ON  THE  HOMOGENIZATION  OF  CLAYS 


By  Horia  I.  Eiie  and  Bogdan  Vernescu 


fn  studying  clays,  an  intensive  research  has  been  done  in  two  directions:  a 
inicroscopic  phenomenological  study,  on  one  hand,  and  a  macroscopic  constitutive 
modeling,  on  the  other.  The  present  research  intends,  by  using  the  homogenization 
techniques,  to  make  a  link  between  the  two  theories. 

In  modeling  the  microscopic  behavior  of  clays  it  is  signiheant  to  consider  the 
presence  of  negative  surface  charges  caused  by  imperfections  of  the  crystal  lattice  of 
the  clay  particles.  The  unbalanced  charges  are  compensated  by  the  accumulation  of 
cations  from  the  solution  surrounding  the  mineral  particles.  A  diffuse  double  layer 
is  created. 

The  problem  is  described  at  the  microscopic  level  by  the  balance  of  momentum 
in  both  parts  of  the  porous  medium  (the  fluid  part  and  the  solid  skeleton) 

—  div<7^  ~  /,  and  —  divo-^  =  /  (1) 

where  the  constitutive  equation  for  the  fluid  phase  involves  an  electric  stress  part: 


' 

’-'ij  ~ 
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with  being  the  stress  tensor  in  the  fluid  in  the  absence  of  the  electric  field. 
We  have  to  adjoin  the  balance  of  mass: 

divu  I) 

the  Maxwell’s  laws: 

div(7/j5)  =  q  in  the  fuid,  d\v{riE)  =  0  in  the  solid 
with  q  the  free  charge  density  and  r;  the  electric  permeability, 

curl/i;  =  0 

and  the  conservation  of  free  charge  in  the  fluid: 

fliv./  =  0 


with  the  free  current  deiisitv: 


J  rjE  'V  qv 


u  being  the  conductivity. 

The  boundary  conditions  are: 

[<7'n]  =  0,  [t;-7i]  =  0,  [u,]  =  ;i/((7^t)r.  (8) 

where  Ur  =  n  -  (u  •  n)n. 

We  have  to  adjoin  to  these,  the  appropriate  boundary  conditions  for  the  electric 
Held  and  free  current: 

[£■,]  =  0,  [}]E-n]-(i\  J-7i  =  0  (9) 

with  the  surface  cliarge  density. 

It  is  important  to  note  that  the  first  jump  condition  in  (8)  reduces  to  the 
double  layer  condition,  that  is  usualy  used  to  describe  the  presence  of  the  electric 
charges  on  the  surface  between  tire  fluid  and  a  rigid  solid. 

The  third  condition  in  (8)  reduces,  in  the  absence  of  the  volume  and  surfaces 
charges,  to  the  slip  condition  used  in  the  case  of  composite.?  with  imperfect  interfaces. 
For  the  case  of  small  viscous  stress  compared  to  the  Maxwell  stress  part  (a°n)r  « 
-  <7’^)7i)r  the  condition  reduces  to  the  C  potential  condition  in  the  Debye-Hiickel 
approximation. 

The  homogenized  equations  display  a  coupled  non-linear  behaviour. 


Schemas  de  volumes  finis  pour  les  problemes 
d’ecoulements  en  milieu  poreux. 

Robert  Eymard^  Thierry  Gallouet^  et  Raphaele  Herbin^ 


Les  problemes  d’ecoulements  en  milieu  poreux  possedent  de  nombreux  aspects 
physiques,  induisant  des  comportements  varies,  de  transports  convectifs  et  dif- 
fusifs  notamment.  Pour  resoudre  numeriquement  ces  problemes,  des  schemas 
numeriques  de  type  ’’Volumes  Finis”  sont  frequemment  employes.  Nous  presen- 
tons  ici  les  principes  generaux  permettant  de  batir  ces  schemas,  fondes,  en 
particulier,  sur  un  principe  de  conservativite.  Nous  donnons  ensuite  quelques 
resultats  de  convergence  des  solutions  approchees  donnees  par  ces  schemEis  pour 
certains  problemes  simples  d’ecoulements  en  milieu  poreux. 
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Boundary  layer  corrections  in  the  poroelastic 
seabed  Lying  under  a  Shallow  Ocean 

Bob  GILBERT 

University  of  Delaware,  Newark,  U.S.A. 

We  consider  the  problem  of  acoustic  wave  propagation  in  a  shallow  ocean,  with  a 
reactive  seabed.  We  shall  model  the  seabed  as  a  poroelastic  structure  using  the  method 
of  homogenization.  The  effect  of  the  ocean  will  be  taken  into  account  by  computing 
boundary  layer  terms.  It  will  be  assumed  that  the  ocean  occupies  the  region 

IR[0,h]  ■■=  {(■^^1,-T2,3:3)  :  0  <  .T3  <  h,  {x^.,X2)  e 

whereas  the  seabed  occupies  iR[o,-oo]  •=  {(^i 1 3:^21  •'1^3)  :  — oc  <  X3  <  0,  (xi,X2)  G  iR  }. 
The  poroelastic  structure  in  the  seabed  will  be  modeled  by  assuming  a  periodic 
microstructure  dependent  on  a  small  parameter  e  which  is  of  the  same  order  of 
magnitude  as  the  pore  size.  In  iR[o,-oo])  shall  refer  to  the  region  as  that 
part  occupied  by  the  solid,  elastic  matrix  and  as  that  part  occupied  by  the  fluid. 
Moreover  we  assume  that  all  the  pores  are  filled  so  that  iR[o,_oo]  —  U  D/f.  We 
consider  the  solid  to  consist  of  an  elastic  porous  material,  which  we  model  in  the  solid 
space-time  region  flse  X  (0,r)  by  using  the  elastic  equilibrium  equation 

and  in  the  fluid  space-time  region  x  (0,T)  using  the  viscous  fluid  equilibrium 
equation 

=  (2) 

^  dt  dxj 

Here  u  :=  (tfi,U2,U3)  is  a  solid  displacement;  whereas  a  (ui,  112,^3)  is  a  flux  (or 
velocity).  The  fluid  stress-tensor  is  given  by 

af;  :--p^^,j+2^(e)e.y(ir).  (3) 

where  is  the  pressure,  /2(e)  a  viscosity  which  depends  on  e,  and  6ij[v  )  a  rate  of  strain 
tensor.  In  addition,  we  assume  that  the  fluid  is  incompressible;  hence 

V  •  =  0  in  (D/e  U  iR[o,/i])  x  (O.T).  (4) 

The  elastic  coeflicients  of  the  solid  matrix  are  taken  to  be  anisotropic.  The  boundary 
layer  corrections  in  the  seabed  due  to  the  ocean  abutting  upon  the  poroelastic  structure 
are  calculated  following  the  scheme  used  in  [1]. 
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THEORICAL  AND  NUMERICAL  DERIVATION  OF  POLYMER 
FLOW  THROUGH  TWO  MEDIA 

O.  GIPOULOUX,  A.  MIKELIC 


After  a  short  review  on  classical  isotherm  viscosity  laws  (Power  law,  Carreau 
law  and  his  variants),  we  propose  the  theorical  derivation  of  the  microscopic 
polymer  flow  in  two  particular  cases  to  the  different  filtration  laws  depending 
on  the  injection  velocity  or  the  Reynolds  number. 

At  first,  we  consider  an  incompressible  bidimensionnal  polymer  flow  through 
a  thin  layer,  by  classical  homogenization,  we  determine  three  filtrations  laws 
wich  are  function  of  the  order  of  the  injection  velocity  in  regard  of  the  thin 
layer’s  depth.  When  the  injection  velocity  is  small  enough  compared  with  the 
thin  layer’s  depth,  we  obtain  an  explicit  relation  between  the  filtration  velocity 
and  the  pressure  gradient,  when  the  injection  velocity  is  great  enough  compared 
with  the  layer  depth,  we  obtain  an  explicit  filtration  law  connecting  the  filtration 
velocity  and  a  certain  power  of  the  pressure  gradient.  Between  this  two  Ccise, 
we  obtain  an  implicit  transition  filtration  law.  Some  numerical  experiments  are 
proposed  to  illustrate  this  results. 

In  a  second  part,  we  consider  the  case  of  a  polymer  flow  through  a  porous 
medium.  After  a  review  of  the  filtration  laws  obtained  by  theorical  homoge¬ 
nization,  we  try  to  propose  some  numerical  algorithms  to  solve  the  homoge¬ 
nized  problems  wich  present  a  coupling  between  the  microscopic  scale  and  the 
macroscopic  scale  due  to  the  two  scale  convergence. 


AN  EXAMINATION  OF  A  STOCHASTIC-NONLOCAL  THEORY  FOR 
MODELING  DISPERSION  IN  SCALE-DEPENDENT  POROUS  MEDIA 


Matthew  M.  Perkins  (Purdue  University,  School  of  Chemical  Engineering) 
John  H.  Cushman  (Purdue  University,  Department  of  Agronomy) 

Robert  A.  Greenkom  (Purdue  University,  Department  of  Chemical  Engineering) 


A  stochastic-nonlocal  theory  used  to  predict  mean  concentration  in  megascopic  systems 
was  investigated.  The  investigation  involved  numerically  solving  the  stochastic-nonlocal 
theory  and  comparing  the  results  to  experimentally  measured  laboratory  data.  The  theory 
uses  the  Darcy  scale  transport  equation  along  with  the  log-hydraulic  conductivity  concept. 
The  velocity,  log-hydraulic  conductivity  and  concentration  are  decomposed  into  their 
mean  and  fluctuation  terms  and  the  resulting  mean  and  mean  removed  equations  are 
solved  by  use  of  Fourier  and  Laplace  transforms.  A  nonlocal  constitutive  flux  is 
incorporated  into  the  theory.  A  FORTRAN  program  was  used  to  numerically  solve  the 
theory.  The  resulting  concentration  profiles  and  second  spatial  moments  were  compared 
to  data  from  a  one-dimensional  tracer  experiment  characterized  by  nonlocal  dispersion. 
The  experiment  involved  a  cylindrical  column  packed  with  glass  beads  in  a  heterogeneous 
fashion.  It  was  found  that  the  simulation  of  the  stochastic-nonlocal  theory  did  not 
adequately  model  the  experimental  results.  The  simulation  concentration  profiles  did  not 
match  the  experimental  concentration  profiles.  Specifically,  the  simulation  significantly 
overpredicted  the  size  of  the  experimental  mixed  zone.  Further,  the  simulation  did  not 
model  the  nonlocal  behavior  well  since  the  results  matched  the  local  advection-dispersion 
equation  better  than  the  nonlocal  experimental  results.  It  was  also  found  that,  while  the 
comparison  of  the  concentration  profiles  was  poor,  a  comparison  of  the  second  spatial 
moments  of  the  concentration  profiles  was  quite  good.  Two  reasons  for  the  poor 
performance  of  the  simulation  are  hypothesized:  1)  the  simulation  solves  for  the  mean 
concentration  while  the  experiment  corresponds  to  a  single  realization  of  concentration;  2) 
the  covariance  function  chosen  for  the  log-fluctuating  conductivity  may  not  actually  model 
this  property  very'  well. 


A  Model  of  Multicomponent  Diffusion 
in  a  Reactive  Acid  Soil 
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Ulrich  HORNUNG  (UniBwM, 
ulrich@informatik.iinibw-muenchen.de) 
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In  structured  soils,  a  considerable  amount  of  soil  water  can  be  immobile. 
As  this  immobile  water  is  confined  to  the  smallest  pore  classes,  it  may  be 
nevertheless  in  contact  with  most  of  the  reactive  surface  areas  of  the  solid 
phase.  A  model  is  proposed  for  the  following  general  case:  i)  The  mobile  soil 
water  is  characterized  by  a  relatively  low  pH  and  high  AP"^  concentrations 
due  to  internal  or  external  inputs  of  acids;  ii)  the  solid  phase  consists  of 
silicate  minerals  that  weather  irreversibly  and  relatively  independent  of  pH 
by  releasing  base  cations,  and  iii)  all  soil  water  is  exposed  to  a  high  partial 
pressure  of  CO2. 

In  the  simplest  case  we  model  this  situation  as  a  one-dimensional  mul¬ 
ticomponent  diffusion  problem  in  which  the  connection  to  the  mobile  phase 
and  the  solid  phase  are  described  by  Dirichlet  and  Neumann  boundary  con¬ 
ditions,  respectively.  Precipitation  of  secondary  minerals  (A1(0H)3)  in  the 
water  column  is  linked  to  decreased  diffusion  coefficients  of  dissolved  ions. 

The  resulting  set  of  coupled  non-linear  partial  differential  equations  has  a 
nonstandard  form,  since  a  reaction-diffusion  system  is  coupled  to  a  condition 
of  electric  neutrality  everywhere  in  space.  Therefore,  special  care  has  to  be 

taken  when  applying  numerical  techniques  for  parabolic  systems. 

Test  calculations  show  that  the  model  exhibits  a  new  mechanism  for 
accumulation  in  acid  soils  that  interferes  in  a  novel  way  with  variations  in 
NOj  in  the  mobile  water. 


TRANSPORT  EQUATIONS  FOR  FLUID  DISPLACEMENTS  IN 
HETEROGENEOUS  POROUS  MEDIA :  THE  M.H.D.  MODEL 


Roland  LENORMAND 

Institut  Frangais  du  Petrole 
BP  311 

F-92  506  Rueil-Malmaison  Cedex 


CPU  time  for  numerical  simulations  can  be  reduced  by  using  improved 
transport  equations  in  an  equivalent  "homogenised"  medium  instead  of 
fine-scale  simulations.  Most  of  simulators  use  the  Koval  non-linear 
equation  for  large  viscosity  ratios  and  homogeneous  media.  We  have 
derived  the  average  equations  to  be  used  when  the  medium  is  highly 
heterogeneous,  and  especially  layered.  Our  method  is  based  on  a  two- 
fluid"  approach  and  the  determination  of  a  matrix  of  relative 
permeabilities  that  accounts  for  anisotropy  and  gravity  effects.  The 
originality  of  the  method,  called  M.H.D.,  is  to  combine  the  numerical 
calculation  of  the  streamtubes  for  unit  viscosity  ratio,  in  any  of  the 
three  main  directions,  to  the  stochastic  calculation  of  the  displacement 
inside  the  streamtubes. 

The  form  of  the  transport  equations,  and  the  values  of  the  parameters 
that  govern  the  equations  are  found  to  depend  strongly  on  the  flow 
regimes.  For  the  limiting  cases  : 

•  dispersive  regime:  a  standard  dispersion  equation  is  valid  with  a 
dispersion  coefficient  "D"  proportional  to  the  correlation  length  and 
permeability  variance  (macrodispersion). 

•  fingering  regime  (due  to  viscosity  or  gravity)  :  the  transport  equation 
reduces  to  the  known  Koval's  equation  with  an  effective  viscosity  ratio 
"M". 

•  for  channeling,  at  unit  mobility  ratio,  it  is  shown  that  the  Koval  non 
linear  equation  is  still  a  good  approximation,  with  a  heterogeneity 
parameter  "H"  instead  of  M. 

A  general  equation  is  then  proposed  for  displacements  in  reservoirs 
with  long-range  correlations  and  viscosity  effects.  The  equation 
combines  the  dispersive  term  (parameter  D)  and  the  non-linear 
convective  term  with  the  product  HM. 

The  advantage  of  the  method  is  that  the  streamtubes  and  the  D  and  H 
factors  are  determined  from  a  monophasic  flow.  It  can  also  be  used  for 
calculating  the  pseudo-permeabilities  for  immiscible  displacements  with 
negligible  capillary  effects. 


CALCULATION  OF  POROUS  MEDIA  EFFECTIVE  PROPERTIES: 
COMPUTATIONAL  PROBLEMS  AND  REQUIRED  UNIT  CELL  FEATURES 

by  P.  FABRIE*,  M.  QUINTARD**,  S.  WITHAKER*** 
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ABSTRACT 


The  calculation  of  effective  transport  properties  starting  with  the  description  of 
the  porous  medium  in  terms  of  a  periodic  unit  ceU  is  particularly  attractive.  It  provides 
bounds  or  even  estimates  that  can  be  very  accurate  and  useful  for  modeling  purposes. 
Besides  the  theoretical  aspects  associated  with  the  derivation  of  the  local  problems  that 
are  needed  to  compute  the  effective  properties,  two  problems  are  of  a  particular 
importance.  The  first  problem  is  the  implementation  of  effective  algorithms,  which  can 
be  a  very  difficult  task  if  convective  effects  are  important,  and  if  the  geometry  is  very 
complicated.  The  second  problem  is  associated  with  the  degree  of  complexity  required 
in  the  unit  cell  geometry  in  order  to  reasonably  represent  actual  characteristics  of  "real" 
porous  media.  In  this  paper,  several  algorithms  are  presented  to  solve  for  local 
problems  associated  with  dispersion-like  processes.  High  accuracy  schemes  can  be 
implemented  by  solving  pseudo-transient  problems  with  operator  splitting.  Details  and 
recommendations  are  provided. 

Resulting  numerical  models  are  used  to  test  for  the  importance  of  the  unit  ceU 
geometry  and  the  type  of  process  under  consideration,  i.e.,  diffusion  or  dispersion.  It  is 
shown  that  highly  dispersive  systems  calls  for  very  complex,  disordered  systems. 
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Upscaling  of  Dispersivity  in  Modeling  of  Solute  Transport: 

Mathematical  Theory  and  Simulations  of  Laboratory  Experiments 
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One  need  in  the  current  theory  of  subsurface  transport  in  porous  media  is  an  im¬ 
proved  understanding  of  upscaled  transport  physics  in  highly  heterogeneous  subsurface 
environments.  This  understanding  should  take  the  form  of  mathematical  models  that 
are  valid  at  multiple  scales.  This  investigation  first  focuses  on  a  Lagrangian  approach, 
similar  to  a  formulation  of  Dagan,  for  calculating  the  dispersion  tensor  coefficients  of 
the  transport  equation.  The  tensor  is  defined  in  terms  of  the  second  spatial  moment 
of  a  plume  with  respect  to  the  plume’s  centroid.  This  form  of  the  dispersion  tensor 
is  derivable  from  the  transport  equation.  The  Lagrangian  approach  characterizes  the 
dispersion  tensor  as  a  stochastic  process,  which  in  turn  leads  to  consideration  of  the 
transport  equation  as  a  stochastic  partial  differential  equation.  The  lecture  discusses 
this  theory  and  numerical  models  using  it  to  simulate  flow  and  transport  in  an  exper¬ 
imental  tank  packed  with  sands  of  varying  hydraulic  conductivities.  Comparisons  to 
experimental  data  are  presented. 


Analysis  of  Finite  Volume  Methods 
and  Application  to  Reservoir  Simulation 

J.-M.  THOMAS  &  D.  TRUJILLO 
Universite  de  Pau,  FRANCE 

A  new  Mixed  Finite  Volume  Method  for  elliptic  problems  on  quadrilateral  meshes 
will  be  presented  and  analyzed.  We  obtain  so  a  conforming  and  conservative  convergent 
scheme,  with  a  9-point  stencil  in  2-D  (27-point  in  3-D).  Generalization  with  a  triangular 
unstructured  mesh  is  then  given.  These  method  is  well  adapted  for  domain  decomposi¬ 
tion.  Some  numerical  examples  in  reservoir  simulation  and  developed  by  Elf-Aqultaine 
Society  will  be  presented. 
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STUDIES  IN  THE  ASYMPTOTIC  REGIMES  OF 
DISPLACEMENTS  IN  POROUS  MEDIA 

Z.  Yang  and  Y.C.  Yortsos 
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University  of  Southern  California,  Los  Angeles,  CA  90089-1211 

ABSTRACT 

We  study  the  asymptotic  behavior  of  miscible  displacements  in  porous  media  of  2-D 
rectilinear  geometry  in  the  two  limits  where  a  permeability-modified  aspect  ratio  Rl  becomes 
large  or  small,  respectively.  In  either  case,  the  problem  reduces  to  the  solution  of  an  integro- 
differential  equation.  We  investigate  the  validity  of  the  asymptotic  description  and  derive 
appropriate  matching  conditions  in  regions  where  it  is  not  valid.  For  the  case  of  unstable 
displacement,  the  linear  stability  characteristics  of  the  equation  as  well  as  the  non-linear 
finger  evolution  are  described.  Comparison  with  the  numerical  solution  of  the  full  problem 
delineates  the  range  of  validity  of  the  asymptotic  regimes. 


On  a  Model  of  Fractured  Porous  Media 


by 

loana-Andreea  ENE*  and  Jeannine  SAINT  JEAN  PAULIN* 


In  this  paper  we  study  a  double  porosity  model  in  a  double  periodicity  media. 
From  a  mechanical  point  of  vue  this  model  represents  a  fractured  porous  media. 
From  a  mathematical  point  of  vue  we  study  a  Neumann  problem  with  double  peri¬ 
odicity.  We  prove  existence  and  unicity  for  such  a  problem  and  using  the  three-scale 
convergence  we  obtain  the  homogenized  equation  and  the  homogenized  coefficients. 
From  a  mechanical  point  of  vue  the  result  we  obtain  is  a  Darcy  law  at  the  macro¬ 
scale. 

At  least  in  the  steady  case  these  result  show  us  that  both  the  double  periodicity 
model,  introduced  by  Th.  Levy  [3]  and  P.  Donato  and  J.  Saint  Jean  Paulin  [2],  and 
the  double  porosity  model,  introduced  by  T.  Arbogast,  J.  Douglas  and  U.  Hornung 

[1],  are  the  same. 
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Non-linear  Darcy’s  Law 
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We  study  the  stationary  Navier-Stokes  system  in  1-dimensionaI  porous  medium. 
Let  V  C  [0, 1]^  be  periodic  in  ei  direction  with  upper  and  lower  boundaries 
r+  and  T~  (dist(r+,r“)  >  0).  For  £  =  1/n  ,  n  €  N  we  define  the  domain 


Let  S  =  U  n  {xi  =  0}  =  y  n  {xi  =  1).  Let  g'  G  Co~(S)  ,  i  =  1,2  be 
such  that  =  J-^g^  0-  For  gi{x)  =  g\x/e)  we  consider  the  problem 

J  -Au= -1- -f  =  0  ,  divu''=0  infl^ 

I  =  0  on  r+"  ,  =  \giei  on 

for  =  eSb  Using  the  asymptotic  expansions  in  the  form 
u"  =  ,  p"  =  \{v^{xi)  +  ep\xuy)  +  ---]  (2) 

we  obtain  the  two-scale  homogenised  problem  - 

f  -A,u“  +  KV>°  +  V,N  =  -gei  ,  div,u°  =  0  in]0,l[xy 
\  =  0  on  r+“  ,  is  1-periodic  in  j/i  ,  fy  u°  =  0  . 

For  \0\  <  yi^/2\/2  we  prove  that  the  problem  (3)  has  the  unique  solu¬ 
tion  u°  e  L2(0,l;  VU(y))  ,  G  L2(]0,l[xy)/R  ,  p°  G  L2{0,1)/R  where 


W{Y)  =  ^  =  0  ,  ^|r+-  =  0  ,  (j)\s  1-periodic  in  j/i}.  More¬ 

over  =  vP{y)  €  W{Y)  ,  y  =  P^{y)  S  L‘^{Y)  and  p'^  is  the  linear  function 
p°{xi)  =  -axi  -f  /5. 

By  computing  boundary  layers  uh  on  E?  ,  i  =  1,2  and  estimating  the  difer- 
ence  we  obtain  the  following  result: 

Theorem  1  Let  Cp{\u°\L^  -f  -|-  |u^,|i,oo}  <  1.  Then 

\eu^  -  [u°(f)  +  u\i{^)  +  . 

jeV-{p°-f  e[pMf)+Pj\(f  )+p^;(f)]}|L2(n^)  <  Ce-i/"  . 

As  the  consequance  we  conclude  that 

/sq-)  •  •  ^2)dx2  fy  u°  weakly  in  L‘^(0, 1) 

£^/sq-)P^(  ■ -^2)^x2 strongly  in  E^(0,1)  , 

where  S^(j/i)  =  e2(yi)  and  S(t)  —  Yf]  {yi  —  t}. 

Finally  we  find  the  relation  between  filtration  velocity  v  —  fy  u°  and  the 
pressure  gradient  ^  in  the  form  of  non-linear  law 


V  =  k{  —  - 


where  k{a)  =  fy  Ui{a)  for  u{a)  E  W(Y)  the  solution  of  the  problem 

J  —  Au(q;) -f  (u(Q)V)u(a)  +  Vp(a:)  =  aci  ,  div?i(Q:)  =  0  in  y  ,  , 
1  u{a)  =  0  on  ,  {u{a),p{a))  is  1-periodic  in  j/i  .  ^  ^ 

For  the  function  k  we  find  the  Taylor’s  expansion  in  neighbourhood  of  0  with 
^(™')(q)  —  Jy  u^^\a)  where  u^”^(a:)  are  the  solutions  of  the  linear  problems 
obtained  by  formal  derivation  of  (7)  with  respect  to  a.  We  prove  that  A:(0)  = 
A:”(0)  =  0  so  that 

— ■  >■> 


moreover  ^'(0)  is  the  standard  Darcy’s  permeability  (as  in  the  linear  law) 
and  A:’”(0)  >  0. 
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by 

D.  Loggia,  N.  Rakotomalala  and  D.  Salin 

Laboratoirc  Fluid.es  Automatique  et  Systemes  TberniK^ues, 

Batiment  502,  Campus  Universitaire,  91405  Orsay  Cedex,  France 

and 

Y.C.  Yortsos 

Department  of  Chemical  Engineering  and  Petroleum  Engineering  Program, 
University  of  Southern  California,  Los  Angeles,  CA  90089-1211,  USA 

ABSTRACT 

When  a  fluid  of  lower  viscosity  displaces  a  fluid  of  higher  viscosity,  the  interface  between 
the  two  fluids  is  unstable.  The  resulting  interface  pattern  is  refereed  to  as  viscous  Angering 
U  When  a  denser  fluid  sits  on  the  top  of  a  lighter  one,  the  situation  is  also  unstable  leading 
to  gravity  fingering  U  When  gravity  and  viscosity  act  on  simultaneously,  they  can  contribute 
either  in  the  same  way  to  the  interface  stability  (  both  stabilizing  or  destabilizing  )  or  in 
opposite.  In  the  latter  case,  two  situations  have  to  be  described.  A  lighter  and  less  viscous 
fluid  is  injected  from  top  to  bottom  against  a  denser  and  more  viscous  fluid  (flow  1):  at  low 
flow  rate,  q  ,  the  gravity  stabilizing  effect  overcomes  the  viscous  destabilizing  one  leading  to 
an  overall  stability.  Above  a  large  enough  flow  rate  qi  viscous  effects  prevail  leading  to  an 
overall  instability.  In  the  reverse  situation  (flow  2)  the  denser  and  more  viscous  is  injected 
from  the  top  against  the  other  one  in  which  case  large  flow  rates  lead  to  an  overall  stability 
due  to  viscous  effects.  Below  a  small  enough  flow  rate  q2,  gravity  effects  prevail  resulting 
to  instability.  In  the  pioneering  approach  of  Hill  ^  the  linear  stability  analysis  (long  waves. 


LW)  predicts  a  similar  threshold  of  instability  for  the  two  reverse  flow  qf^  =  =  q^,  in 

reasonable  agreement  with  the  experiment  with  small  gravity  and  viscosity  contrasts. 

qh^kg— - —  (1) 

If  the  LW  ap roach  is  relevant  for  immiscible  fluids  interfaces  where  the  interface  thickness 
(b)  can  always  be  considered  as  small  compared  to  the  wave  length  (A)  of  the  the  perturbation 
involved  in  the  stability  analysis,  this  no  longer  true  for  miscible  fluids  because  of  longitudinal 
dispersion  which  spreads  the  interface  as  times  goes  on  (6  ~  y/i).  Therefore  the  details  of 
the  front  concentration  profile  can  be  felt  by  short  waves  {  X  b)  leading  to  new  predictions 
For  miscible  fluids  a  non-monotonic  viscosity  or  a  gravity- viscosity  profiles  can  leads  to 
SW  instabilty  whereas  LW  theory  predicts  stability.  These  new  instability  branches  in  the 
dispersion  relation  leads  to  qij^q^- 


ft  = 

In  the  poster,  we  discuss  the  state  of  the  art,  present  a  derivation  of  this  new  thresholds  of 
instabilty  and  compute  accurately  the  dispersion  relation,  i.e.  growth  rate  of  the  disturbance 
versus  wave- vector.  For  the  numerical  computation,  we  address  the  case  of  real  fluids  for  the 
concentration  viscosity  and  density  relationships  to  get  results  suitable  for  comparison  with 
experiments 


l_Mar  ^9  l_A/m 


(2) 
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On  considere  un  domaine  fluide  obtenu  en  retirant  d’un  ouvert  fl  des  parties  solides  eYg 
par  periodicite.  Pour  e  fixe,  on  considere  le  systeme  de  Stokes 


Trouver  ,p^)  €  x  L‘^(Qc)  / 

-Au-'  +  Vf  =  / 
div(u- )  =  0 
u"  =  0 


et  de  Navier-Stokes 


Trouver  {v'  ,q^)  G  77^ (nt)"  x  / 

—e'^Av‘  +  (u‘.V)u-  +  =  / 

div([;‘ )  =  0 
V'  —  Q 


dans 
dans  fie- 
sur  50^ 


dans  Qe 
dans  Of 
sur 


ou  /  est  une  fonction  appartenant  a  et  0  <  7  <  3/2.  On  definit  ensuite 

et  par 


— =  e-’  dans 
div(w7’^)  0  dans  fie 

—  0  sur  dQc  \  fl 


div(y.0=)  ^  dans  fl,. 

y',i,£  _  g  sur  50^  \  fl 


avec  A',  matrice  de permeabilite  de  composante  j  On  considere  le  probleme 

suivant,  qui  va  permettre  de  definir  la  vitesse  de  Darcy  ou  vitesse  homogeneisee  dans  fl.. 

(  Trouver  G  R'(fl)  / 

I  div(A'(/- Vp°))  =0  dansfl 

[  A'(/- =  0  sm  dQ 

On  definit  ensuite  la  vitesse  de  Darcy  par  v  =  K{f  -  Vp°).  Nous  rappelons  que  si  A  est  un 
ouvert  de  R",  alors  on  note  \a.  la  fonction  caracteristique  de  I'ensemble  A.  On  a  alors 

Theoreme  0.1  Soit  {u'.p")  la  solution  du  probleme  de  Stokes  (1).  On  note  R-  =  {r  G 
fl  /  dist{x,dQ)  >  y/s]  et  u-  le  prolongement  par  0  sur  D  de  u' .  On  suppose  que  f  £ 
On  a  alors  les  resultats  de  convergence  suivants: 


XB'nn,- {p' -P°}  — >0  dans  L-{n) 


ou  p^  est  fixe  par  la  relation  -  Po)  =  0 


HO 


Nous  avons  egalement  un  resultat  similaire  pour  le  probleme  de  Navier-Stokes.  Cependant, 
la  convergence  de  la  pression  n’est  acquise  qu’avec  une  restriction  supplementaire  sur  7. 
Ainsi,  pour  des  valeurs  de  7  proches  de  3/2  nous  allons  etre  oblige  de  faire  appel  a  la 
construction  d’un  operateur  de  relevement  pour  la  pression. 

Theoreme  0.2  Soit  (v'  ,q-)  la  solution  du  probleme  de  Navier-Stokes  (2).  On  suppose  que 
f  €  et  que  7  verifie  0  <  7  <  |  si  n  =  3  et  0  <  7  <  |  sz  u  =  2.  En  notant  v,  le 

prolongement  de  par  0  sur  fi,  on  a  les  resultats  de  convergence  suivants: 

-^0  dans 

J  'J  ' 

Si  7  est  tel  que  0  <  7  <  |  si  n  =  3  e<  0  <  7  <  1  si  n  =  2,  alors  on  a 


XB’nn,- iq' -  >  ■ — J- 0  dans  L^{Q) 


XB'nn, •  {?' -P°}  ^  0  dans  L~(0) 


ou  q^  est  fixe  par  la  relation  {q~  —  Po)  —  0. 

Si  7  est  tel  que  |  <  7  <  |  si  n  =  3  et  1  <  7  <  |  si  n  =  2,  alors  il  existe  r  6]l,  2[  et  un 
relevement  q^  de  q^  appurtenant  a  Lq{Q)  tel  que 


q'  —  p°  —  e  ^  — >  0  dans  ^5(0) 

J 

q'  —  0  dans  Ln{0) 

£-(■0  “ 
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The  uranium  mining  activities  and  remediation  of  the  deposits  in  the  chalk  formation  of 
Straz  of  the  Northern  Bohemia  is  the  main  subject  of  the  contribution  proposal.  Nowadays, 
the  mining  activity  is  being  suppressed  and  the  problems  concerning  the  remediation  of  the 
deposits  in  Straz  appear  to  be  very  urgent.  One  of  the  most  important  aspect  which  is  face 
is  the  protection  of  the  detached  roof  stratum  -  the  reservoir  of  the  drinking-water  m  this 


is  necessary,  most  of  aU.  to  consider  much  larger  region  than  it  is  in  the  current  models 
and  to  compute  the  transport  of  contaminants  applying  three-dimensional  models^  With 
respect  to  the  large  size  of  the  problems,  necessity  of  non-standard  use  of  the  FEM  and 
particularities  of  uranium  mining  in  the  considered  region,  it  is  possible  to  exploit  the  alrea  y 
known  experience  only  partially.  Thus  the  whole  process  of  creating  the  new  models  is  stiU 
open  The  solution  of  this  problem  is  divided  into  two  parts:  In  the  first  part  -  the  solution 
of  stationary  porous  media  fiow  of  the  technological  solution.  It  is  assumed  that  the  chemical 
non-homogenity  of  the  liquid  does  not  affect  its  porous  media  flow.  In  the  second  part,  the 
transport  of  chemical  substances  in  the  technological  solution  and  the  chemical  changes  in 
the  rock  are  studied.  Both  particular  problems  are  solved  by  several  computational  methods. 

The  solution  of  underground  water  flow  problem  in  the  real  conditions  must  reflect  com¬ 
plex  .Geological  structure  of  sedimented  minerals.  Layers  of  the  stratified  rocks  with  substan- 
tiaUy°different  physical  properties  must  be  modelled  using  the  appropriate  discretization  o 
the  geological  region.  These  geological  characteristics  can  be  correspondingly  described  by 
the  mixed-hybrid  finite  element  method  using  trilateral  prismatic  elements  with  vertical  faces 

and  generally  nonparaUel  bases.  rr^i.  •  i 

Let  U  be  a  bounded  domain  in  with  a  Lipschitz  continuous  boundary  dQ.  The  potential 

fluid  flow  in  a  saturated  porous  media  can  be  modeUed  by  the  velocity  u  given  by  Darcy  s  law 
u  _  where  p  is  the  piezometric  potential  (fluid  pressure)  and  A  is  symmetric 

and  uniformly’ positive  definite  second  rank  tensor  of  hydrauhc  permeability  of  the  porous 
medium.  Consider  also  the  continuity  equation  for  the  incompressible  flow  V  •  u  -  q,  where  q 
represents  the  density  of  potential  sources  in  the  medium.  The  boundary  conditions  are  given 
by  p  =  PD  on  u-n  =  UjY  on  OQn,  where  dQ  =  SUcUSOiv  are  such  that  -  0 

and  n  is  the  outward  normal  vector  defined  (almost  everywhere)  on  the  boundary  dH. 

Assume  from  now  that  the  domain  D  is  polyhedron  and  is  subdivided  in  a  collection  ol 
subdomains  £h,  such  that  every  subdomain  is  a  trilateral  prism.  We  denote  the  coUection  of 

‘This  work  was  supported  by  GA  CR  under  grant  201/93/0067. 


faces  of  subdomains  e  e  €h  which  are  not  adjacent  to  the  boundary  dflo  by  T/i  =  - 

SfijD,  where  h  is  the  discretization  parameter. 

The  velocity  function  u  will  be  approximated  with  the  vector  functions  linear  on  every 
element  e  e  Sh  and  the  pressure  function  p  will  be  approximated  with  the  element-wise  con¬ 
stant  function  and  the  trace  of  the  pressure  on  the  structure  of  faces  Th  will  be  approximated 
by  the  face-wise  constant  function. 

Hybrid-mixed  formulation  of  the  porous  media  flow  problem  leads  to  the  system  of  linear 
equations: 


In  the  second  part,  the  transport  of  chemical  substances  in  the  technological  solution 
and  the  chemical  change  in  the  rock  are  studied.  The  hybrid-mixed  model  of  the  transport 
of  chemical  substances  dissolved  in  technological  solution  is  described  by  Pick’s  law  = 
-DV'c,  where  'i  is  the  flow  vector  of  1-th  substance  and  ‘c  denotes  its  concentration,  D  is 
the  diffusivity-dispersivity  tensor  calculated  by  the  formula  [D]ij  =  DmSij+aT-\n\.8ij  +  {aL- 
^  where  Dm  denotes  coeficient  of  molecular  diffusion  and  oj  respectively  ai  denotes 
coefficient  of  transversal  dispersivity  respectively  coefficient  of  longitudinal  dispersivity.  The 
equation  of  the  materials  balance  for  1-th  substance  is  ^  -f  u  •  V^c  +  'V  -‘i  +  ^cQ  —‘r{c,t)  = 
,  where  first  member  defines  storage  of  1-th  substance,  second  member  defines  convection 
and  third  expresses  the  influence  diffusion  and  dispersion.  The  member  ^cQ~  determines 
the  influence  exhausting  of  1-th  substance  and  the  member  ‘c*Q'^  determines  the  influence 
indentation  of  1-th  substance  about  concentration  'c*.  Last  member  on  the  left  side  equation 
defines  the  complicated  chemical  reaction  of  components  of  the  technological  solution  and 
component  of  the  rock  by  that  1-th  substance  either  rises  or  extinguishes. 

This  member  is  analysed  using  reaction  kinetics  and  decomposed  on  the  part 
where  coefficient  ^R~{t)  determines  a  measure  of  a  capacity  of  the  environment  bring  into 
reaction  with  1-th  substance  and  on  the  part  ‘R'^{t).  which  determines  a  measure  of  a  capacity 
of  the  environment  generated  1-th  substance. 

Balance  equations  for  aU  substances  are  discretized  in  the  time  with  the  step  At  =  and 
space  on  the  equal  decomposition  £h  of  the  domain  fi  as  in  the  porous  media  flow  problem. 
The  finite- dimensional  approximation  of  the  transport  chemical  substances  problem  using 
the  function  spaces  described  above  leads  to  the  following  system  of  equations: 

/  D  B  C\  /‘in\  (  '51  \ 

‘Hn  j  j  =  j  ,  ('-ffn)('cn)  =  ('52,n), 

where  superscripts  I  —  denote  substances  in  the  solution  and  superscripts  I  =  L\  -j- 

1, ...,  L  denote  substances  in  the  rock. 

We  perform  the  numerical  experiments  with  the  problem  of  transport  of  chemical  sub¬ 
stances  is  testing  at  present  time  using  a  set  of  basis  problems. 
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La  prise  de  moyenne  a  rendu  possible  I'obtention  de  I'equation  de  transport  macroscopique 
en  milieux  poreux,  Dans  le  cas  du  probleme  de  dispersion,  I'equation  macroscopique  est  similaire 
a  I'equation  microscopique  si  les  variables  locales  sont  remplacees  par  les  variables  moyennes  et  si 
le  tenseur  effectif  de  dispersion  est  introduit  a  la  place  du  coefficient  de  diffusion  moleculaire.  Le 
tenseur  peut  etre  calcule  si  on  connait  le  champ  des  vitesses  a  I'interieur  de  la  structure  poreuse  et 
si  on  est  en  mesure  de  resoudre  le  probleme  de  fermeture  associe  a  la  prise  de  moyenne.  Le 
probleme  de  fermeture,  est  un  probleme  classique  de  convection-diffusion.  Cependant,  la 
geometrie  du  milieu  est  une  grande  inconnue  et  le  choix  que  Ton  fait  pour  aborder  numeriquement 
la  resolution  du  probleme  de  dispersion  en  milieux  poreux  peut  etre  determinant  pour  la  reussite  de 
la  methode.  La  question  qui  se  pose  est  de  savoir  si  Ton  peut  trouver  des  arrangements 
geometriques  simples  qui  puissent  representer  raisonnablement  un  milieu  reel.  Le  choix  c'est  porte 
ici  sur  des  arrangements  periodiques  bidimensionnels  de  cylindres  carres.  Pour  certains  auteurs,  si 
la  periodicite  spatiale  autorise  une  resolution  numerique,- voire  analytique,  aisee  du  probleme  de 
dispersion,  elle  engendre  necessairement  des  comportements  singuliers  du  tenseur.  Une  etude 
numerique  complete  de  revolution,  en  fonction  du  nombre  de  Peclet  des  composantes 
longitudinale  et  transversales  du  tenseur  de  dispersion  est  proposee  dans  des  arrangements 
d'abord  simples  (en  ligne,  en  quinconce)  puis  plus  complexes  (disposition  aleatoire  de  la  phase 
solide  dans  la  cellule  unite).  On  s'interessera  plus  particulierement  a  I'influence  de  la  direction  0  de 
I'ecoulement  moyen  par  rapport  aux  axes  de  la  cellule  unite. 

Pour  des  milieux  dits  "ordonnes",  le  coefficient  de  dispersion  longitudinal  presente,  pour  0=0°  et 
45°,  une  variation  proche  de  Pep2 ,  assez  differente  de  celle  obtenue  pour  d'autres  valeurs  de  0. 
Pour  le  coefficient  de  dispersion  transversal,  les  resultats  sont  fortement  dependants  de 
I'arrangement  considere. 


Pour  les  20  milieux  "desordonnees"  etudies,  les  resultats  sont  quasi  independants  de  Tangle  0  et 
les  variations  moyennes  sont  en  Pep^>6  et  Pep^’"^  pour  le  coefficient  respectivement  longitudinal  et 
transversal. 

Les  resultats  pour  les  milieux  "ordonnes"  sont  confrontes  a  ceux  issus  de  la  theorie  approchee  de 
Koch  et  Brady  dont  on  pent  dire  qu'elle  est  uniquement  valable  lorsque  la  porosite  du  milieu  tend 
vers  Tunite.  Cette  conclusion  est  confortee  quand  on  aborde  la  comparaison  des  resultats  obtenus 
ici  avec  ceux  issus  des  etudes  de  Salles  et  al.  puis  Edwards  et  Eidsath  et  al.  pour  des  arrangements 
de  meme  type. 

Dans  le  but  d'elargir  la  comparaison  a  des  resultats  experimentaux,  peu  presents  dans  la  litterature, 

un  dispositif  experimental  est  mis  au  point.  L'ensemble  est  compose  de  deux  dispositifs  :  le 

premier  a  pour  objet  la  fabrication  de  milieux  poreux  modeles  par  photo  polymerisation  laser.  Le 

principe  de  fabrication  repose  sur  une  technique  appelee  stereophotolithographie  laser  developpee 

par  le  GRAPE  (ENSIC  NANCY)  permettant  de  realiser  des  objets  de  tres  bonne  definition 

geometrique.  Les  milieux  ainsi  obtenus  sont  formes  d'un  arrangement  periodique  de  cylindres  de 

section  circulaire  de  5  mm  de  hauteur  et  de  0,5  mm  de  diametre  disposes  entre  deux  plaques  de 

verre.  Un  milieu  comporte  en  moyenne  4000  obstacles  repartis  sur  une  largeur  egale  a  25  mm 

pour  un  miheu  de  longueur  egale  a  130  mm.  Les  experiences  de  dispersion  reposent  sur  la  mesure 

de  revolution  de  la  concentration  d’un  traceur  colore  introduit  dans  Tecoulement  au  sein  du 

milieu.  Cette  mesure  est  realisee  par  un  ensemble  video,  camera  et  analyseur  d'images.  Les  profils 

temporels  de  concentration,  prise  en  valeur  moyenne  sur  une  periode  spatiale  du  milieu,  analyses 

par  une  methode  d'ajustement  sur  le  modele  theorique,  permettent  d'estimer  la  variation  du 

coefficient  de  dispersion  longitudinal  en  fonction  du  nombre  de  Peclet.  L'interet  majeur  de 

l'ensemble  experimental  est  que  les  resultats  obtenus  peuvent  etre  compares  de  fagon  effective  aux 

resultats  numeriques  donnes  pour  des  geometries  identiques.  Les  premiers  resultats  concement  un 

arrangement  en  hgne  etudie  pour  deux  directions  de  Tecoulement  moyen  :  0°  et  26°.  L’exposant  m 
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de  la  loi  limite  pour  le  coefficient  de  dispersion  longitudinal  =  A  Pep"'  presente  bien  une 

forte  sensibilite  a  la  direction  de  Tecoulement,  conformement  aux  predictions  numeriques.  Le 
decalage  sur  les  valeurs  effectives  du  coefficient  (inferieures  a  0°,  superietnes  a  26°,  comparees 
aux  valeurs  obtenues  numeriquement)  semble  devoir  etre  attribue  a  la  forme  de  Tobstacle.  Un 
milieu  desordonne  est  etudie  pour  une  direction  du  vecteur  vitesse  moyenne.  La  cellule  unite  du 
milieu  etudie  experimentalement  est  identique  a  Tune  des  20  geometries  aleatoires  etudiees 
numeriquement.  Dans  ce  cas  Taccord  sur  les  valeurs  effectives  du  coefficient  de  dispersion 
longitudinal  est  excellent  cependant  que  Ton  obtient  m=l,2,  valeur  generalement  admise  pour  les 
structures  poreuses  "reelles",  sensiblement  plus  faible  que  la  valeur  predite  numeriquement  (1,38). 
Les  trois  experiences  realisees  demontrent  Taptitude  du  dispositif  pour  etudier  qualitativement  et 
quantitativement  le  processus  de  dispersion.  L'outil  mis  au  point,  meme  s'il  reste  encore 
perfectible,  couple  a  un  code  de  calcul  Liable  doit  permettre  une  meilleure  comprehension  des 
differents  mecanismes  d'un  processus  particulierement  complexe. 
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Ce  travail  a  pour  but  de  determiner  I'ordre  de  variation  des  termes  non  lineaires  de  la  loi 
de  Darcy,  c'est-a-dire  comment  varie  la  perte  de  charge  en  fonction  du  vecteur  vitesse 
moyenne.  La  plage  de  nombre  de  Reynolds  particulaire  consideree  est  choisie  de  fa9on  a  ce  que 
I'ecoulement  se  trouve  dans  le  regime  laminaire  inertiel  [1].  Differents  types  d'arrangements 
periodiques  de  cylindres  carres  ont  ete  retenus  :  arrangements  disposes  en  ligne,  en  quinconce, 
en  chicane  et  de  fagon  aleatoire  pour  une  porosite  de  0,6  environ.  L'etude  s'attache  a  fake 
varier  le  nombre  de  Reynolds  particulaire  et  la  dkection  du  vecteur  vitesse  moyenne  de 
I'ecoulement  par  rapport  aux  axes  du  maillage. 

La  methode  des  volumes  finis  est  adoptee  pour  resoudre  I’equation  de  Navier-Stokes  en 
regime  permanent  avec  des  conditions  aux  limites  periodiques  sur  la  frontiere  du  domaine  de 
resolution,  c'est-a-dke  une  cellule-unite  du  milieu  poreux  periodique.  Deux  types  de  schemas 
numeriques  sont  utilises  pour  la  discretisation  des  termes  convectifs  de  I'equation  ;  un  schema  a 
trois  points  et  un  schema  quadratique  a  quake  points  afin  de  kaiter  la  diffusion  numerique.  Les 
termes  diffusifs  sont  discretises  a  I'aide  d'un  schema  des  differences  cenkees.  La  resolution 
utilise  I'algorithme  SIMPLE  (Semi-Implicite  Method  for  Pressure-Linked  Equations)  developpe 
par  Patankar  et  Spalding  [2]. 


Les  resultats  montrent  que  pour  les  maillages  et  critere  de  convergence  retenus,  les 
valeurs  calculees  sent  sensiblement  equivalentes  pour  les  differents  schemas,  ceci  pour  des 
nombres  de  Reynolds  particulaires  inferieurs  a  400  environ. 

Les  resultats  cherchent  plus  particulierement  a  decrire  au-dela  du  regime  de  Stokes  les 
deviations  obtenues  par  rapport  a  la  loi  de  Darcy  a  savoir  la  correction  a  apporter  a  la  perte  de 
charge  lineaire  ainsi  que  la  perte  d'isotropie  observee  pour  les  arrangements  periodiques 
reguliers  possedant  deux  axes  de  symetrie  geometriques  perpendiculaires.  Les  variations  de  la 
perte  de  charge  nous  amenent  a  conclure  a  I'existence  de  deux  types  de  correction  non  lineaire  a 
la  loi  de  Darcy.  Pour  une  plage  intermediaire  de  nombre  de  Reynolds  particulaire,  la  correction 
non  lineaire  est  d'un  ordre  superieur  a  une  variation  quadratique  relativement  a  la  vitesse 
moyenne  <u>,  proportionnelle  a  <u>n  ou  n  est  superieur  a  2.  Une  valeur  de  I'exposant  n=3 
semble  correcte  sur  une  plage  de  nombre  de  Reynolds  particulaires  moderes  pour  des  milieux 
isotropes  lorsque  les  termes  d'inertie  sont  importants  mais  ne  sont  pas  preponderants  devant  les 
termes  visqueux.  Cette  constatation  est  en  accord  avec  les  resultats  theoriques  prevus  par 
Wodie  et  Levy  [3],  Mei  et  Auriault  [4]  et  Moura-Neto  et  al  [5]. 
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Gianni  Baldinf*^ 


I  the  experiments  I  performed  in  the  lahorator)'  of  p.p.A., 

Trieste)  (Itjtly),  in  order  to  explain  what  happens  during  the  water  ^rcolation.  I 
made  experinsents  with  water  temperature  of  4''C,  since  1  want- to  the 

effects  of  extraction  process.  As  qualitative  results,  I  have  that  the  flow  pte  9, 
oorrespohdlhg  to  a  fixed  Ap,  is  decreasing  \rith  time  and  that  the  asyniptQtiC  valuc 
of  the  flow  h  hot  necessarily  an  increasing  functloD  of  Ap  (1].  I  Supp^  there  are 
two  main  pfoc^ses  related  with  'this  behaviour:  first  of  all  the  finest  dompopent-  of 
the  solid  petioles  is  removed  from  the  coarse  skeleton  and  tr^SpOrted  W  the 
water  flow,  U  ^that  it  tends  to  form  a  compact  layer  whose  hydraulic  _^nductivity 
is  consider^ly  smaller;  on  Ih^  other  hand  1  have  a  (partMy  Jdrrovcrslble) 

■  d^ormation’of  the  skeleton  of  the  porous  medium  induced  by  the  flow.  The  first 
one  is  investlgUed  in  [2],  where  the  proof  of  the  weD^j>osedness  ^  the  ^problem.,  of 
the  unicity  df  the  solution  and  of  the  continuous  dependence  upon  the  data  and  the 
coefficients  Is  j^ven;  moreover,  it  is  shown  in  (2J  that  there  arc  sufficient  conditions 
for  a  finite  tioic  extinction  of  the  transport  phenomenon  and  that  its  occurrence  is 
indeod  necessaty  in  order  that  the  asymptotic  discharge  can  exhibit  a  nonmohotone 
dependence  upon  the  external  pressure.  Concerning  the  model  of  the  deformation  of 
the  pkeleton'i  if  d{t)  denotes  the  thickness  reduction  of  the  coffee  bed  at  time  t,  I 
can  relate  it  to  the  water  flow  9;  for  instance  by  assuming  that.  ^  ^  f}q{0k"^)  +  , 

'  - 


Matheni^ticaJ  niodeis  for  freezing  in  porous  media 

Fedehico  Talamuoci 
/ 

If  v>'e  observe  a  column  of  porous  soil  saturated  vith  water  and  subject^  to  ^  freezing 
temperature  at  tbe  top  (the  lower  side.'on  the  contrary,  is  maintained  at  6  positive  .tempera^ 
ture),  we  will  notice  fundamentally  two  facts:  a  frtezing  fmnt  njotdng  downw^d  (tonfining 
temperature  permitting)  end  a  water  migration.' irorsii  the  bottom  towards  higher  parts  of 
the  soil.  The  mass  transfer  is  due  to  capillarity  phenomena,  or  rather  to  cfyoiuction. 

locate  essentially  tw*o  parts  of  the  soil:  an  unfrozen  region,  with  a  porous  inatri>c 
and  liquid  water,  and  a  i^ozen  part,  v^diere  ice  fills  the  porous  space.  |  ' 

■  Nevertheless,  what  can  happen  in  proximity  to  the  separation  frofit  is;  the  formation 
of  pure  ice  layers,  called  lences,  that  are  fed  by-  the  migration  flux  of  water  toT^ds  the 
freezed  part.  If  we  consider  that  thedower  base  of  the  soil  is  at  rest,  we  wdll  Observe  a 
lifting  of  the  upper  surface,  due  to  the  formation  of  lences:  such  a  process  is  known  as 
fro.fthoave,  ;  i  '  ; 

On  the  other  hand,  if  the  freezing  process  is  too  fa?t,  or  an  ovcrU/Wnifn  (a 

load  that  the  column  must  support)  pievents  the  soil  fiom  separating  somewhere’,  wo  will 
simply  observe  a  downward  movement  of  the  freezing  front  {frott  penetraUon).'  A  froti 
heave  can  be  observed  even  in  this  case,  due  to  the  increase  of  the  mass  of  Water  (the  soil 
IS  maintained  saturated):  Instead  of  a  macroscopic  lens,  a  change  in  porosity  occur, 
after  the  passage  of  the  frost  front.  ■  1 

'  If  we  assume  that  the  separation  fregion  between  the  two  main  p^ts  df  the  soil  is  a 
surface  {freo^^mg  front) ^  we  are  dealing  with  a  class  of  models  which  are  more  treatable 
from  a  mathematical  point  of  view*.  In  this  case,  the  lifting  of  the  soil  is  knowp  a?  primary 
frost  heave.  ; 

Nevertheless,  the  experimental  observation  of  the  process  makes  lis  rertiark  the  exis¬ 
tence  of  an  intermediate  region,  situated  between  the  frozen  and  the  unfrdzeu  part  and 
called  frozen  fringe^  where  ice  and  ■v.’ater  coe.xist  in  the  porous  space.'  Models  including 
this  transition  zone  describe  the  secondary  frost  heave.  1 

In  any  case,  the  starting  point  for 'writing  the  equations  of  any  model  is  the  conserva¬ 
tion  of  heat  and  mass.  Obviously,  we  need  further  information  about  the  evolution  of  the 
process.  ' 

.  Three  questions  arc  of  crucial  importance:  i 

—  which  is  the  .composition  of  the  frozen  fringe,  that  is,  which  are  the  inost  relevant 

interfacial  effects  in  the  mixture  iCe- water- soil?  .  : 

—  which  is  the  law  regulating  the  water  flux  in -the  frozen  fringed  \ 

—  are  temperature  and  pressure  iudojpcndc-nt  in  the  same  region? 


Models  proposed  in  literature  differ  one  from  llic  other  just  on  the  giound  of  the  kind 
of  answer  given  to  the  previous  questions;  the  thermodynamical  and  transfer  process  in 
the  fro7/en  fringe  are  still  something  obscure  and  there  is  not  a  convergence  of  opinions 
among  modelers  about  the  govenring  laws. 

The  analysis  that  we  are  doing  deals  with  a  quasl-stcady  approach  to  the  problem, 
assuming  that  the  temperature  is  linear  with  respect  to  the  spatial  coordinate  (we  are 
considering  one  dimensional  models). 

As  regards  primary  fro hcav^,  we  found  out  conditions  on  the  (given)  confining  tern* 
peralui^B  such  that  we  are  able  to  distinguish  between  lens  formation  Kud  frost  pentiraUen 
Farthermore,  a  proof  of  an  existence  theorem  has  been  concluded  in  the  general  case  (not 
linear  temperature). 

The  investigation  of  secondary  fivsi  heavt  is  started  through  the  analysis  of  a  particular 
model  proposed  recently  by  Y.  Nakeno:  he  assumes  that  temperatme  arid  water  PJ'essurc 
are  independent  in  the  frozen  fringe  and  the  hydraulic  flux  is  regulated  by  a  combinado^,, 
of  the  thermal  and  pressure  gradients.  j  | 

Moreover,  the  volumetric  v.-ater  content  in  the  frozen  frinyt  is  assumed  to  be  a  known  j 
function  of  the  temperature.  Ice  is  at  rest  with  respect  to  the  porous  matrix  (in  opposition 
to  that  models  whero  ice  moves  as  one  rigid  body)  and  always  at  the  atmospheric  pressme 
(instead  of  a  time  and  space  dependence  of  the  icc  pressure,  responsible  in  that  case  for 
the  soil  separation  as  soon  as  lens  starts  to  form). 

In  particular,  a  criterion  which  allows  us  to  discriminate  between  the  two  possibil¬ 
ities  wDl  be  analyzed  and  the  existence  of  the  solution  In  a  general  case  (v/ithout  any 
simplification)  will  be  investigated. 
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1  Introduction 

This  work  originates  from  an  attempt  to  model  from  mathematical  and  numerical  points 
of  view,  the  leaching  of  ions  from  cement  matrices  used  for  waste  solidification  and  waste 
storage.  The  aim  of  this  model  is  to  study  the  leaching  of  ions  on  a  large  scale  of  time 
which  can’t  be  done  by  experimental  measurements.  Here,  we  present  a  simple  model  for 
leaching  of  several  ionic  chemical  species  accounting  for  phase  change,  diffusion  and  chemical 
reactions. 

Some  numerical  results  are  discussed  in  the  case  of  one  chemical  species  (lime)  and  two 
chemical  species  (lime  and  lead).  We  present  a  modification  of  the  original  model  to  fit  with 
precipitation  phenomena. 

2  Mathematical  model  for  diffusion  and  dissolution 


We  consider  the  process  of  diffusion  and  dissolution  on  the  macroscopic  scale.  Interstitial 
solution  of  cement  is  assumed  to  be  diluted  and  at  rest.  The  quantities  studied  are  Cjt  the 
molar  concentration  of  ionic  chemical  species  k  in  liquid  phase  and  Sp  the  molar  concen¬ 
tration  of  chemical  species  p  in  solid  phase.  We  denote  l{p)  the  set  of  indexes  of  species  k 
yielding  from  dissolution  of  species  p  in  solid  phase.  The  porosity  0  <  ^  <  1,  is  assumed 
to  be  constant. 

The  domain  of  the  study  is  a  cylinder  fi,  filled  with  hydrated  cement,  during  a  period 
oftime[0,T]. 

fi  =  |(x,y,  r)  G  R^,  +  <  1,  0  <  z  <  l|  (2-1) 

Using  the  conservation  principle,  we  consider  the  variation  of  the  relative  number  of 
moles  of  species  k  filling  an  arbitrary  subdomain  U  of  Q.  This  variation  is  given  by  the  sum 
of  the  diffusive  flux  through  the  boundary  and  the  internal  production/destruction  rate. 

According  to  Pick’s  law,  the  outgoing  diffusive  flux  of  ionic  species  k,  qk  through  the 
surface  dV  (the  boundary  of  subdomain  V)  expresses  as: 

qk=DkVCk-n  1  <  ^  <  (2-2) 

where  the  constant  Dk  stands  for  the  diffusive  coefficient  and  n  for  the  outgoing  normal 
to  the  volume  V. 


3  PHASE  CHANGE  RATES 


We  denote  by  RkiCi,  ..,Cn)  the  molar  rate  of  production  of  species  k  due  to  chemical 
reactions  in  liquid  phase,  and  for  I  <  k  <  L  we  denote  by  Fk{Ci,  ..,Cn>  Sp)  the  molar 
production  rate  of  species  k  yielding  from  dissolution  of  species  p  in  solid  phase. 

Hence  we  have; 

^  Jy  6Ck  dxdydz  =  fgy  6Dk  VCk  ■  nda+ 

Jy  dRk(Ci, Cn)  dxdydzP  (2-3) 

fy(l  -  0)Fk(Ci, Cn,  Sp)  dxdydz,  l<k<N. 

By  using  Green’s  formula,  and  since  equations  2-3  are  valid  in  every  smooth  subdomain 
V ,  we  have: 

edtCk  -  div{0DkVCk)  =  9RkiCu  ..,Cn)  +  (1  -  d)Fk{Ci,  ..,CN,Sp)  l<k<N.  (2-4) 

It  remains  to  specify  the  dissolution  kinetic  equations  which,  because  of  mass  balance 
between  solid  and  liquid  phases  read: 

=  -  ^  F,(Ci,..,Cn,5p)  1<p<M.  (2-5) 

keiip) 

One  of  the  major  assumptions  is  that  chemical  reactions  in  liquid  phase  are  at  equilib¬ 
rium,  thus  Ejfe  =  0  for  1  <  A:  <  A’;  Moreover,  the  model  is  supposed  to  be  one-dimensional. 

For  0  <  T  be  given,  the  equations  2-4  and  2-5  become: 


OdtCk  -  d,{dDkd,Ck) 
dtSp 


(1  -  d)Fk{Ci,  ..,Cn,Sp)1  <  k  <  N,  in(0,l)x(0,T)  (2-6) 

-  Y,  Fk{Cu..,CN,Sp)l<p<M,  in(0,l)x(0,r)(2-7) 

*£/(?) 


3  Phase  change  rates 


The  main  phenomenon  involved  is  dissolution  of  solid  phase  in  open  system  (see  [3]). 
There  is  no  general  expression  in  literature  to  express  the  dissolution  rate  in  such  a  case.  In 
this  model  we  consider  it  is  proportional  to  the  difference  between  the  concentration  of  the 
species  in  solution  Ck  a.nd  a  reference  concentration  C^ ,  provided  that  solid  phase  is  still 
present: 


r  Fk  =  0iov  L  +  1  <k  <  N-  /o 

\  FkiCi, ..,  Cn,  Sp)  =  osffn+(5p)[Q:  -  Ck]+  for  1  <  A:  <  I  with  0  <  o.  ^  ^ 

Cl  stands  for  thermodynamical  equilibrium  concentrations  of  ionic  species  k.  It  is  deter¬ 
mined  by  a  system  of  equations  including  equilibrium  relationships  for  species  in  solution, 
the  electrical  equilibrium  relation  and  the  equilibrium  relationships  for  dissolution.  The  last 
ones  are  replaced  by  an  equality  between  Ck  and  C^.  when  there  is  no  more  of  the  related 
solid  phase. 

Existence  and  uniqueness  of  solutions  to  the  mathematical  problem  for  a  single  species 
in  solid  phase  have  been  proved  in  a  previous  article  (see  [1]).  A  numerical  method  based  on 
a  marching  technique  to  approximate  numerically  the  mathematical  problem  for  diffusion 
and  dissolution  was  also  proposed  (see  [2]).  This  method  was  first  tried  in  the  case  of  a 
single  species  and  extended  to  the  case  of  two  species. 

In  the  case  of  two  species,  computations  were  done  with  lead  and  lime.  The  results 
showed  that  we  reached  values  of  pH  where  usely  a  phenomenon  of  precipitation  appears. 
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This  is  due  to  the  fact  that  we  took  into  account  two  dissolved  species  for  lead.  The 
precipitation  rate  is  not  written  our  model. 

To  consider  the  possibility  of  precipitation,  the  model  is  changed  as  follows: 

^k{Ci,  ..,Cn,  Sp)  —  aisgn^ {Sp)[C).  —  Ck]^  —  o:2[Cf;  —  Ck]  ,  1  <  ^  <  h.  (3"9) 

Numerical  simulations  with  a  given  shows  that  Ck  ’’follows”  the  equilibrium  concen¬ 
tration  and  that  precipitations  can  occur  with  this  model.  This  expression  is  interesting 
because  it  represents  both  dissolution  and  precipitation  and  because  it  may  be  use  for  en- 
vironemental  flows  too  (see  [4]).  Results  got  with  this  model  about  the  example  lime  and 
lead  will  be  discussed. 
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ABSTRACT 

Lattice  gas  has  been  proved  useful  to  simulate  flow  as  an  alternative  to  Navier- Stokes 
finite  difference  method.  We  will  demonstrate  the  relevance  of  such  a  simulation  tool  in  the 
case  of  flow  of  miscible  fluids  in  a  Hele-Shaw  cell.  This  is  a  particular  and  simple  way  to 
approach  the  properties  of  viscous  flow  inside  porous  media  (Darcy’s  law  holds).  Using  a 
recent  and  efficient  lattice  gas  model  we  simulate  the  mixing  at  the  interface  of  two  viscous 
miscible  fluids  in  such  a  flow,  and  adress  the  question  of  the  resulting  instabilities. 

To  simulate  the  flow  we  use  the  lattice  gas  BGK  model  proposed  by  Qian  et  al.  [1] 
that  describes  the  fluid  dynamics  in  2D.  This  model  is  an  alternative  to  conventional  lattice 
Boltzmann  models  ([2],  [3])  which  themselves  evolved  from  the  lattice  gas  automaton  [4]. 
The  basic  variable  in  the  model  are  V,(r,t)  which  denotes  the  mean  occupation  number 
of  particles  and  Aj-(r,  t)  which  represent  the  relative  concentration  of  the  two  fluids  in  the 
direction  i  at  the  node  r  at  time  t.  These  number  are  updated  according  to  the  Boltzmann 
equation 


W(f  +  Q,/  +  1)  =  Ni{rA)  +  ARiV,(r,Q  -  ivr(r- t))  (1) 

and  a  convection- diffusion  equation 

Ai(f  +  qW  +  1)  -  A,-(r,Q  +  AD(Ai(r,t)  -  A^^(f,i))  (2) 

where  is  the  relaxation  parameter  that  determines  the  kinematic  viscosity  v  and  A^ 
the  one  who  determines  the  molecular  diffusivity  Dm  and  c,-  denotes  the  velocit}"  vector 


connecting  neighboring  nodes  on  the  lattice.  In  2D  we  take  a  square  lattice  with  eight 
directions  of  propagation  (including  the  diagonals)  plus  one  at  rest  (model  D2Q9).  In  3D  a 
cubic  lattice  is  used  with  fifteen  or  nineteen  directions  {D3Q15  or  D3Q19  models). 

The  algorithm  involves  three  steps:  given  an  initial  uniform  distribution  of  Ni  and  the 
chosen  initial  amount  of  color  A,-,  the  equilibrium  distributions  N^'^{r,t)  and  A®’(T, t)  are 
calculated;  then  a  collision  step  (right  hand  side  of  equation  1,2)  is  performed;  finally,  the 
occupation  numbers  are  propagated  to  the  neighboring  sites.  In  case  the  last  site  belong  to 
a  solid  boundary  (top  or  bottom  walls),  they  are  propagated  back  in  the  opposite  direction 
to  enforce  a  no-slip  boundary  condition  at  the  walls.  Along  the  flow  direction  x,  periodic 
boundary  conditions  are  applied  for  the  Ni  as  for  the  A,-.  Finally,  to  simulate  an  imposed 
flow,  a  fixed  small  amount  e  of  the  mean  occupation  number  on  each  node  of  the  first  column 
of  lattice  is  transferred  from  the  upstream  to  the  downstream  direction  at  each  time  step. 
The  simulations  are  performed  on  a  256  x  32  lattice  for  the  2D  studies  and  on  a  256  x  32  x  8 
lattice  for  the  3D  studies. 

On  2D  we  perform  simulation  with  a  real  2D  flow,  i.e.  in  a  channel,  using  both  the 
original  Navier-Stokes  and  convection-diffusion  equations  (1,  2).  On  the  other  hand  we  test 
the  Rayleigh- Taylor  instability,  where  the  flow  is  induced  by  gravity.  Finally  we  build  a  3D 
Hele-Shaw  cell  with  our  lattice.  Our  results  are  compared  to  classical  numerical  simulations 
[5]  and  to  an  Hele-Shaw  experiment  performed  in  our  laboratory.  The  effects  of  the  viscosity 
ratio  are  presented,  for  the  stable  case  and  the  instable  one.  The  application  to  more  complex 
systems,  as  real  porous  media,  will  be  now  considered. 
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Nonlinear  diffusion  in  an  inhomogeneous 

aquifer 


M.  Guedda,  D.  Hilhorst  and  M.  Peletier 


We  consider  the  Cauchy  problem 

f  p{x)ut  ■=  (u(  1  -  u)Ua:)x  111  jR  X  jR+ 

[  u{x,  0)  =  Uo{x)  X  £  R. 

We  sui>pose  that  the  functions  p  and  uq  satisfy  the  hypotheses 

p  £C\R)rM^{R)r\L^{R),  p>0 


uo  £  0  <  uo  <  l,no(l  —  '^o)  has  compact  support. 

The  related  equation  p{x)ut  =■  Au™'  has  been  studied  by  Kamin  and 
Rosenau  and  Kamin  and  Kersner.  In  particular  they  prove  that  in  the  case 
that  the  space  dimension  n  >  3  the  support  of  u  may  become  unbounded 
in  finite  time  if  p  decreases  fast  enough  as  \x\  tends  to  infinity.  Our  main 
purpose  is  to  prove  a  similar  property  for  Problem  (P). 

After  having  shown  how  Problem  (P)  modehzes  the  interface  between 
fresh  and  salt  groundwater  in  an  inhomogeneous  aquifer,  we  sketch  the  proofs 
of  existence  and  uniqueness  of  a  weak  solution  u  of  Problem  (P). 

We  then  give  a  result  about  the  limiting  behavior  of  n  as  t  tends  to 
infinity. 

Theorem  A  (Large-time  behaviour)  Let  u  be  the  solution  of  Problem 
(P).  Then 

_  j  p{x)uo{x)dx 

u[t)  — )■  u  =  ^  . -  as  t  — >  oo, 

/  p(x)dx 

J  R 

uniformhj  on  compact  subsets  of  R. 


as  t  — >  oo, 


A  solution  u  is  said  to  exliibit  finite  time  blow-up  (FTB)  if  there  exists  a 
time  T  such  that  suppu(t)(l  —  u(t))  is  unbounded  from  above  for  all  times 
t  >  T.  For  the  formulation  of  the  next  result  we  also  need  an  auxiliary- 
density  function  cr  defined  by 

a{x)  =  min{p(^)  :  0  <  ^  <  x}, 

the  reason  being  that  the  function  a  is  monotonic  while  p  need  not  be. 

Our  main  result  is  the  foUo-wing. 

Theorem  B  Let  u  be  the  solution  of  Problem  (P).  Then  the  following  im¬ 
plications  hold: 

fiOO 

1.  If  /  xp{x)  dx  <  oo,  then  finite  time  blow-up  occurs: 

Jo 

poo 

2,  If  x(r{x)  dx  =  oo^  then  there  is  no  finite  time  blow-up. 

Jo 


W.  Anker,  D.  Graillot,  J.  Bourgois  et  M.  Zelfani : 

Spatialisations  quantitative  et  qualitative  en  milieu  poreux 
Application  a  I’estimation  des  risques  dans  la  depollution  des  sols 

La  spatialisation  dans  les  methodes  quantitatives 

Les  methodes  quantitatives  sont  des  methodes  largement  utilisees  pour  decrire  les  phenomenes 
de  transfert  en  milieu  poreux  souterrain.  La  resolution  du  probleme  de  transfert  fait  appel : 

(i)  :  a  une  equation  aux  derivees  partielles  et  a  la  connaissance  des  conditions  limites  et  initiates, 

(ii)  :  a  la  dimension  spatiale  du  domaine  pour  une  discretisation  de  I’espace  a  considerer. 

Si  ces  methodes  ne  posent  pas  ou  peu  de  probleme  pour  des  milieux  poreux  homogenes,  il  en 
va  autrement  pour  les  milieux  poreux  heterogenes.  Pour  ces  demiers  cas,  un  maillage  adapte  de 
fagon  optimale  et  resultant  d’une  methode  de  discretisation  aux  differences  fmies  avec  mailles 
variables  ou  aux  elements  finis  est  necessaire  de  telle  sorte  que  chaque  maille  coixesponde  a  une 
partie  du  domaine  relativement  homogene  (nous  exclurons  les  milieux  heterogenes  fissures  : 
milieu  «poreux»  a  grande  echelle).  L’optimisation  du  maillage  discretisant  le  milieu  considere 
souleve  cependant  un  probleme  particulier  :  il  s’agit  de  determiner  le  nombre  et  la  forme  des 
elements  du  maillage  en  function  des  connaissances  du  milieu  (mesures  disponibles,  localisa¬ 
tion  des  heterogeneites). 

Deux  approches  differentes  peuvent  resoudre  ce  type  de  probleme  : 

-  les  modeles  deterministes  qui  prennent  bien  en  compte  les  heterogeneites  mais  qui  demandent 
des  donnees  souvent  inconnues  dans  la  pratique. 

-  les  modeles  stochastiques  qui  foumissent  des  resultats  en  termes  de  probabilite  puisque  les 
donnees  foumies  par  la  geostatistique  le  sont  egalement. 

La  spatialisation  dans  les  methodes  gualitatives 

Les  methodes  qualitatives  de  resolution  ne  font  pas  appel  a  des  donnees  uniquement  numeri- 
ques,  mais  aussi  a  des  donnees  descriptives  qualitatives  ou  semi-quantitatives  (eloigne,  rappro- 
che,  immediat  par  exemple)  et  foumissent  egalement  des  resultats  descriptifs.  Parmi  ces 
methodes,  nous  avons  distingue  les  systemes  a  bases  de  connaissances  et  les  systemes  de  simu¬ 
lation  qualitative. 

La  simulation  qualitative  s’interesse  a  la  simulation  des  phenomenes  physiques  par  I’interme- 
diaire  de  raisonnements  sur  leur  comportement  et  sur  leur  changement  caracteristique.  Ces  mo¬ 
deles  s’appuient  sur  des  regies  logiques  qui  relient  differents  parametres  :  «modeles  de 
parametres».  Ces  outils  servent  souvent  a  decrire  des  systemes  qui  se  foimalisent  difficilement 
par  les  methodes  conventionnelles  (ebullition  d’un  liquide  par  exemple).  De  nombreux  pheno¬ 
menes  physiques  font  appel  a  une  dimension  spatiale  qui  a  ete  introduite  par  certains  auteurs 
dans  la  simulation  qualitative,  mais  une  telle  description  se  limite  pour  I’instant  a  des  systemes 
mecaniques  tres  simples.  Si  la  simulation  qualitative  est  appliquee  au  deplacement  d’une  parti- 
cule  dans  un  domaine  discretise,  le  nombre  d’etats  qualitatifs  croit  de  maniere  exponentielle  tres 
rapidement  [86  etats  qualitatifs  pour  4  mailles  !].  Pour  simuler  un  ecoulement,  il  faudrait  de  plus 
tenir  compte  de  I’heterogeneite  du  domaine,  de  la  notion  des  phenomenes  continus,  des  niveaux 
piezometriques,  des  phenomenes  de  transfert  et  de  la  dispersion  dans  le  cas  de  pollution.  Ce  type 
de  methodes  qualitatives,  bien  que  prenant  en  compte  la  dimension  spatiale  pour  des  cas  sim- 


pies,  ne  sont  pas  encore  en  mesure  de  simuler  des  phenomenes  aussi  complexes  que  I’ecoule- 
ment  et  le  transfert  de  polluants  dans  un  domaine  poreux. 

Dans  les  systemes  a  bases  de  connaissances,  la  dimension  spatiale  a  ete  introduite  principale- 
ment  pour  la  construction  de  topologies  et  la  localisation  d’objets  dans  une  topologie.  Les  don- 
nees  et  les  mecanismes  d’inference  sur  les  proprietes  spatiales  peuvent  etre  introduits  dans  les 
systemes  a  bases  de  connaissances  pour  construire  des  regies  de  logique  (par  exemple)  qui  ser- 
vent  de  support  au  raisonnement  spatial. 

Application  a  1’ estimation  des  risques  des  pollutions  dans  la  decontamination  des  sols 

Les  bases  de  connaissances  spatiales  peuvent  etre  utilises  comme  bases  de  connaissances  geo- 
graphiques.  En  realite,  la  complexite  des  phenomenes  de  pollution  met  en  evidence  la  comple- 
mentarite  et  la  necessite  d’associer  les  methodes  de  spatialisation  quantitatives  avec  les 
methodes  qualitatives.  Les  methodes  quantitatives  sont  davantage  destinees  a  simuler  la  propa¬ 
gation  des  phenomenes  physico-chimiques,  les  methodes  qualitatives  sont  davantage  destinees 
a  raisonner  sur  les  resultats  quantitatifs  pour  une  aide  a  la  decision  (mise  en  place  d’un  plan  d’in- 
tervention).  La  contribution  de  chacune  de  ces  deux  categories  de  methodes  est  plus  ou  moins 
importante  selon  le  cas.  Elle  est  en  particulier  guidee  par  I’existence,  la  nature  et  la  distribution 
des  informations  spatiales  (localisation  des  polluants  de  differente  nature,  reperage  des  zones 
vulnerables). 

La  methodologie  proposee  consiste  a  mettre  en  place  un  systeme  d’analyse  spatiale  permettant 
d’evaluer  la  situation  de  la  pollution  en  terme  quantitatif  (scenarios).  A  partir  des  resultats  quan¬ 
titatifs  obtenus,  une  base  de  connaissances  permet  une  evaluation  en  terme  de  risques  ou  d’im- 
pacts  de  la  pollution  (aide  a  la  decision). 


Water  retention  characteristics  for  fractal  soil 
Nigel  Bird 

Silsoe  Research  Institute,  Wrest  Park  .oe,  Beford,  U.K.  MK45  4HS 

For  soil  structure  modelled  as  a  fractal  -  ith  the  pore  space  as  it’s  complement,  it  is  possible 
.’erive  the  cumulative  pore  volume  as  a  ranction  of  pore  scale.  From  this  powerlaw  water 
iC’.cfition  models  may  be  constructed.  However,  pore  connectivity  effects  neglected  in  the 
above  approach  are  likely  to  make  the  relation  between  pore  volume  and  water  retention 
characteristics  more  complex,  and  the  feasibility  of  inferring  fractal  structure  and  fractal 
dimensions  from  experimental  water  retention  data  must  then  be  questioned. 

Pore  connectivity  effects  relate  to  the  accessibility  of  pore  water  to  the  air-water  interface 
within  the  soil  (resulting  in  delayed  drainage),  and  the  accessibility  of  the  pore  water  to  a 
drainage  route  through  the  soil  (resulting  in  entrapped  water). 


In  a  conventional  non-tractal  problem,  the  air-interface  accessibilit}'  issue  diminishes  as 
drainap  progresses  beyond  air  breakthrough.  In  contrast  for  a  fractal  structure,  the 
accessibility  to  an  air  interface  remains  a  problem  at  all  scales  and  thus  at  all  stages  of  the 
drainage  process. 

As  desaturation  of  a  fractal  structure  progresses,  the  probability  that  the  remaining  water- 
filled  porosity  percolates  falls.  If  the  connectivity  of  the  water  filled  porosity  is  disrupted  at 
an  early  stage  of  the  drainage  process,  then  a  substantial  quantity  of  entrapped  water  mav 
result,  and  further  drainage  will  occur  for  only  a  subset  of  the  remaining  water  filled  porosity 
, hydraulically  connected  to  the  water  sink.  This  ignores  the  hydraulic  connections  made  by 
thin  films,  and  the  water  filled  porosity  existing  below  the  fractal  regime,  but  in  practical 
terms  insufficient  time  may  be  allowed  for  these  to  operate. 

These  two  effects  combined  may  result  in  a  nontrivial  dependence  of  water  retention  data 
upon  the  pore  volume  distribution.  Hence  the  interpretation  of  powerlaw  exponents  derived 
from  such  data  in  terms  of  fractal  dimensions  should  be  made  with  caution. 
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This  work  presents  a  progress  report  towards  building  physically  and  structurally  realistic 
models  of  flow  through  porous  media,  based  on  the  local  change  of  scale  method.  The  local 
change  of  scale  technique  leads,  from  first  principles,  to  a  spatially  smoothed  equation  (e.g., 
Darcy's  law)  and  to  a  closure  problem,  set  at  the  microscopic  scale,  over  a  local  Representative 
Elementary  Volume  (R.E.V.)  of  Darcy's  type  (Vp(ro)). 

Vpfro)  must  be  larger  than  the  micro-geometrical  R.E.V.  (Vp(^o))  --  (minimal  volume  of 
length-scale  ^Q  describing  the  geometrical  variability  over  some  length-scale  L ;  ']/q<<^). 

Development  of  a  valid  geometrical  R.E.V.,  VpC'i/o).  coupled  with  its  relationship  to  the 
local  Darcy's  R.E.V.  Vp(ro)  will  allow  explicit  calculation  of  the  permeability  K  as  an  implicit 
function  of  the  micro-geomeUy. 

In  accordance  with  the  existence  of  a  local  closure  form,  the  three-dimensional  micro¬ 
structure  defines  K.  Nonetheless,  realistically  characterizing  the  three-dimensional  micro¬ 
geometry  is  complex,  current  observations  being  of  porosity  exposed  on  slices  intersecting  the 
medium. 

Development  of  a  valid  geometrical  R.E.V.,  Vp('i.o),  and,  a-fortiori,  usefulness  of  the  local 
change  of  scale  method  appears  to  be  dependent  on  one's  ability  to  remove  several  theoretical 
problems  upstream  from  the  local  change  of  scale  model  in  itself: 

-  complete  characterization  of  the  micro-geometry  in  terms  of  N  two-dimensional 
measurable  parameters  n  —  Micro-geometry  =  f({R}i^] «), 

-  reduction  of  the  set  to  a  minimal  subset  relevant  to  the  involved 

physical  property,  namely  K, 

-  development  of  an  "ad-hoc”  procedure  solely  constrained  by  the  subset  (P*},=i;s-> 
allowing  generation  of  synthetic  three-dimensional  porous  media  acting  as  an  interface  between 
the  image-analysis-based  bi-dimensional  knowledge  one  has  of  the  micro- structure  and  the  three- 
dimensional  numerical  codes  based  on  the  local  change  of  scale  technique. 

Herein,  one  intends,  in  two-dimensions,  to  come  up  with  partial  answers  to  the  afore¬ 
mentioned  upstream  problems,  capitalizing  on  a  more  empirical  approach  due  to  Ehrlich  et  al. 
[Ehrlich  et  al.,  1991a  ;  McCreesh  et  al.,  1991].  The  relevance  of  the  plan  to  the  volume  has  been 
empirically  based  by  Ehrlich  et  al.  who  suggest  that  a  finite  set  of  parameters  namely, 

pore  types  and  their  relative  proportions  (obtained  by  deconvolving  the  two-dimensional  micro- 
geometry)  are  adequate  to  characterize  the  three-dimensional  micro-structure: 

N 

Micro-geometry  =  f(^  (Pore  Type‘ ,  Relative  Proportion’ ) ) 

i=l 

Each  pore-type  is  arranged  in  continuous  and  independent  circuits  at  the  origin  of  strong 
correlations  between  pore-types  and  throat  sizes.  Such  relation  between  pore-types  and  throat- 
sizes  has  been  quantified  by  combining  two-dimensional  image-analysis-based  data,  pore  types, 
and  three-dimensional  petrophysical  data,  mercury  drainage.  In  sandstones,  the  tendency  of  pore 
types  to  fill  in  mutually  exclusive  ranges  of  pressure  intervals  would  be  the  consequence  of  the 
proclivity  of  pores  of  like  type  to  be  mutually  adjacent,  connected  by  similarly-sized  throats,  hence 
a  heterogeneous  fabric  that  results  in  a  series  of  quasi-parallel  flow  circuits. 

Using  an  "interpretative"  implicit  model  of  percolation  type  [Ehrlich  et  al.,  1991b],  the 
physically  relevant  geometrical  information  is  reduced  to  a  simpler  sub-set,  {P’}i=ij^={di,NPi)i^i 


to  address  the  relation  k  -  Micro -geometry  under  the  form  k  =  g{di,NP;};^jj^,  where  Npj  is  the 
number  of  pores  of  the  i*  pore  type  per  |im^  and  d^  is  the  diameter  of  throat  of  the  i*^  pore  type  per 
|im. 


Herein,  one  means  to  develop  a  similar,  but  devoid  of  any  constitutive  assumption 
methodology.  The  two-dimensional  micro-structure  is  fully-characterized  by  a  measure  of  type 
{(^}i,,i_3={F(Vx,Vy),(p(Vx,Vy),Ep}  where  F(Vx,Vy)  and  (pCv^.Vy)  are  the  Fourier  modulus  and  phase  of 
the  micro-geometry  and  £p  the  associated  probability  density  function.  The  minimum  subset 
included  in  {F(Vx,Vy),(p(Vx,Vy),£p}  and  held  to  describe  accurately  the  coupling  K  -  Micro¬ 
geometry  by  means  of  the  change  of  scale  method  is  =  {|F(v)p  ,8p}.  |F(v)p  allows 

quantitative  assessment  of  the  natiu-e  and  length-scales  of  the  stmctural  components  of  the  micro¬ 
structure  i.e.,  of  the  minimum  size  of  Yppyp)  wherein  the  micro-geometry  is  stationary  and  ergodic 
(locally  homogeneous). 

A  set  of  synthetic  media  are  generated  from  a  single  set  {|F(v)p;£p}  derived  jfrom  any 
medium  by  convolution  of  these  parameters  with  a  gaussian  field.  Adequacy  of  a  "measure"  of 
type  {|F(v)p  ,£p}  as  a  measure  of  the  two-dimensional  micro-geometry  is  tested  from  an  empirical 
(but  practically  validated)  angle  by  linking  qualitatively  the  stmctural  parameters  {P}i=i  ^  =  {Pore 
typei ,  Relative  proportion*  of  Ehrlich  et  al.  to  characteristics  of  |F(v)P  [Anguy  et  al.{in  press  ; 
Prince  et  al.,  in  press].  This'is  done  by  showing  that  filtered-and-inverse-Fourier-transformed 
images  might  reveal  a  hierarchy  of  micro- stmctural  components  preferentially  containing  distinct 
pore  types  -  throat  sizes  relationships.  The  spatial  arrangement  of  pores  of  like  type,  as  predicted 
by  Ehrlich  et  al.  would  be  the  same  as  that  one  displayed  on  the  filtered  images.  The  amount  and 
nature  of  the  geometrical  information  carried  by  {|F(v)p  ;£p}  are  quantitatively  assessed  by 
comparing  the  Pore-Types  of  real  and  synthetic  media. 


The  local  change  of  scale  method  is  not  absolutely  restricted  to  periodic  media.  However, 
development  of  a  "numerically-tractable"  closure  form  requires  a  weak  periodic  condition 
imposed  on  the  boundaries  of  VQ(ro).  That  is,  K  is  not  calculated  using  the  local  geometrical 
R.E.V.  Vp('i/o),  but  a  local  periodic  Darcy's  R.E.V.,  Vp(ro),  derived  from  Vp(vo)  [Aguy,  1993].  It  is 
commonly  assumed  that  K  derived  this  way  is  a  close  approximation  of  the  permeability  tensor  of 
natural  porous  media,  provided  that  Vp(ro)  is  large  enough.  However,  the  random  character  of 
natural  porous  media  requires  assessment  of  the  ratio  "n"  between  the  size  'i/q  of  Vp(vo)  and  the  size 
ro=nx'i/o  of  Vp(ro).  Because  natural  porous  media  are  not  totally  periodic,  Vp(ro)  must  be  larger  than 
Vp('i/o)  in  order  for  random  effects  to  be  correctly  accounted  for.  On  account  of  practical 
constraints,  this  problem  has  been  so  far  investigated  in  two  dimensions.  A  P*  rough  estimate  of 
"n"  is  derived  using  synthetic  isotropic  media  by  calculating  K  for  unit  cells  of  increasing  sizes 
and  common  centre.  Results  achieved  this  day  verify  that  the  Darcy's  R.E.V.,  Vp(ro),  is  larger  than 
the  geometrical  R.E.V.,  Vp('!/o):  ro=3()x^o  [^'^nguy,  1993].  These  results  are  also  encouraging  in  that 
they  suggest  that  the  realizations  of  synthetic  media  constrained  by  the  same  set  (|F(v)p;8p} 
converge  at  about  the  same  values  of  K.  In  this  respect,  these  results  can  be  taken  as  a  necessary’ 
condition  in  view  of  basing  that  K  is  solely  constrained  by  {|F(v)p;£p}. 
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A  three  dimensional  mathematical  model,  for  flow  and  solute 
or  energy  transport  in  saturated  porous  media  is  presented. 
This  model  allows  to  take  into  account  the  study  of  flow  in  a 
confined  aquifer  or  in  an  unconfined  aquifer;  density  and 
viscosity  of  the  fluid  might  be  considered  as  linearly 
dependant  on  transport  (VOSS  [5] )  ;  during  mass  transport,  the 
miscible  polluant  can  be  chemically  degraded  or  absorbed; 
finally  this  modeling  can  be  realized  in  an  steady  or  unsteady 
way . 

This  model  is  based  on  physical  equations  which  one  are 
described  with  accuracy  in  Bear  [1],  [2]  .  In  this  study,  they 
are  written  in  order  to  consider  particular  conditions  like 
pumping,  infiltration  as  boundary  conditions  of  flux 
(Jacob  [3]).  This  method  allows  to  take  into  account,  for 
example  some  specifications  like  a  pumping  in  upper  layer  of 
the  aquifer.  A  such  problem  induces  to  solve  a  system  with  two 
coupled  partial  differential  equations.  Each  of  them  is  solved 
by  a  finite  element  method.  For  reasons  described  by  MOSE  [4], 
a  two  degree  interpolation  is  .  required  to  get  a  good 
approximation  of  the  velocity,  thus  the  finite  element  used  is 
a  hexahedron  with  20  nodes,  it  permits  to  have  the  previous 
conditions  on  degree  interpolation.  The  finite  element  method 
requires  to  solve  at  each  time  step  some  linear  systems.  In 
order  to  get  a  fast  and  accurate  computation,  a  preconditioned 
conjugate  gradient  method  is  used;  the  same  process  is  used  to 
solve  the  non  symmetric  system  produced  by  solving  transport 
equation . 


When  the  density  or  the  viscosity  depends  on  concentration  or 
temperature,  a  coupled  system  of  two  equations  has  to  be 
solved.  For  this,  an  iterative  method  on  the  approximation  of 
the  pressure  P  is  developed. 

A  description  of  the  model  including  the  physical  equations 
and  the  numerical  schemes  is  presented.  An  example  of  salt 
intrusion  into  coastal  aquifer  is  analysed,  according  to  a 
study  made  with  the  Australian  National  University  of  Canberra. 
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Title:  Use  of  the  Fourier J.aplace  transformation  and  of  diagrammatical  methods  to 
in  heterogeneous  reservoirs. 

Abstract 


interpret  pumping  tests 


Progresses  in  the  power  of  computers  and  in  reservoir  characterization  allow  the  simulation  of 
prcssure-transients  in  complex  reservoirs  generated  with  stochastic  tools.  Classically,  the  interpretation  of 
these  ti'ansicnts  gives  useful  informations  about  the  reservoir  structure  :  a  major  goal  is  to  interpret  these 
tiansients  in  a  stochastic  context.  To  adress  this  problem,  an  important  question  is  to  relate  the  pre.ssure 
variations  at  time  t  to  the  local  permeability  map.  Using  the  Fourier  iranformation  for  spatial  coordinates,  as 
well  as  a  Laplace  transformation  for  time  allows  to  tackle  this  task.  This  formalism  is  used  in  conjuction 
with  a  peiturbaiion  scries  expansion  in  powers  of  the  permeability  fluctuations  to  get  an  explicit  solution. 

ihc  N-ih  order  term  of  this  series  involves  the  hydrodynamic  interaction  between  N  permeability 
heterogeneities. 


The  goal  of  the  work  is  to  ensemble  average  the  pressure  over  all  the  random  permeability  field 
realizations  to  derive  an  equation  obeied  by  the  mean  pressure.  To  adress  this  problem,  Fcyimian  graplis  arc 
introduced  which  allows  a  graphical  interpretation  of  the  perturbation  series.  It  is  shown  that  the  .sum  of  the 
so  Called  one-pailicle  irreducible  graphs  gives  the  kernel  of  a  linear  iiitcgro-differcntial  equation  obeyed  by 
the  ensemble  averaged  pressure,  All  the  information  about  the  heterogeneities  is  thus  contained  in  this 
renormalized  kernel.  We  recover  thus  directly  results  assumed  by  Indclman  and  we  get- a  systematic  and 
direct  scheme  to  evaluate  successive  approximations  of  the  kernel. 

This  equation  iiscll  is  the  starting  point  of  useful  asymptotic  rS’^uUs  and  approximah  ons .  In  particular 
jt  is  shown  that  the  interpretation  of  pumping  tests  yields  the  steady-state  equivalent  permeability  after  a 
sufficiently  long  time  foi  an  infinite  reservoir,  a  rather  intuitive  result.  This  means  that  the  homogcncizaiiou 
process  performed  by  pumping  tests  is  well  understood. 

Jt  is  also  shown  that  the  kernel  of  the  mean  equation  i.s  closely  related  to  the  time  variations  of  some 
spatial  moments  of  the  pressure  field.  These  moments  arc  important  as  they  quantify  the  size  of  the 
in\estigaled  volume  at  any  lime  and  play  a  key  role  for  simplified  pumping  tests  simulation  methods.  This 
equation  is  also  a  good  starting  point  to  study  fractal-like  media  having  coiTclation.s  at  all  Icngthscalc,  In  that 
vase,  wc  can  expect  lime  and  scale  dependant  behaviour  of  the  apjjarcnl  parameters  known  under  the  name 
of  anomalous  diffusion.  Understanding  of  these  phenomena  will  be  of  great  theoretical  and  practical  interest. 
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Abstract :  This  paper  is  devoted  to  the  modelling  of  compressible  flow  through  porous  media. 
In  contrast  with  incompressible  flows,  compressible  flows  through  porous  media  have  been 
little  studied  in  recent  years.  However,  in  addition  to  more  traditional  domain,  new  applications 
of  porous  media,  such  as  the  external  insulation  of  spacecraft  [1],  require  a  proper  modelling  of 
compressible  flows  in  porous  media.  Naturally,  it  exists  a  large  variety  of  compressible  flow 
depending  on  the  Mach  number  and  the  Re>molds  number,  the  coupling  effect  with  thermal 
effects,  etc.  Here,  we  limit  our  scope  to  isothermal  low  Reynolds  and  Mach  numbers  transient 
flows.  Traditionally,  this  type  of  flow  at  the  macroscopic  level  is  described  by  the  classical 
Darcy's  law  combined  with  a  mass  balance  that  includes  the  transient  term.  This  model  is  called 
the  "classic  model"  in  the  follovvdng. 

The  aim  of  this  paper  is  to  check  numerically  the  relevance  of  this  classic  model.  To  this 
end,  we  consider  the  flow  of  an  ideal  gas  within  a  two  dimensional  model  porous  medium 
consisting  of  an  arrangement  of  square  cylinders.  The  flow  is  due  to  the  imposed  pressure 
decrease  at  the  outlet  of  the  fluid  domain.  At  the  microscopic  level,  the  flow  is  computed  by 
solving  the  full  Navier-Stokes  equations  in  two  dimensions.  We  use  the  PISO  (Pressure- 
Implicit  with  Splitting  Operators)  method,  [2].  This  method  is  well  adapted  to  solving  variable 
density  low  Mach  flows.  Special  attention  is  given  to  the  outlet  boundary  conditions  by  using 
the  NSCBC's  method  (Navier  Stokes  Characteristic  Boundary  Conditions,  [3]).  The  results  are 
spatially  averaged  to  obtain  macroscopic  data.  At  the  macroscopic  level,  the  classic  model  is 
solved  by  a  Galerkin  finite  element  method.  The  comparison  between  the  macroscopic  data, 
obtained  on  the  one  hand  by  spatially  averaging  the  microscopic  results  and  on  the  other  hand 
by  solving  the  problem  directly  at  the  macroscopic  level  is  performed  for  various 
depressurization  conditions.  We  exhibit  situations  for  which  a  good  agreement  is  found 
between  the  two  series  of  data  and  situations  for  which  discrepancies  are  observed.  These 
various  behaviours  are  discussed  in  terms  of  the  various  tirrie  scales  controlling  the  flow. 

1  A.Bouhouch,  M.Prat  &  S.Bories.  Transient  compressible  flow  and  thermal 
transfer  within  a  heterogeneous  porous  system.  Journal  of  Thermophysics  and  Heat 
transfer  in  press. 

2  R..I.  Issa,  1986,  Solution  of  the  implictly  discretised  fluid  flow  equations  by  operator¬ 
splitting.  Journal  of  Comp.  Phys.,  Vol.62,  No.l,  pp.40-65 

3  T.J.Poinsot  &  S.K.Lele,  1992,  Boundary  conditions  for  direct  solutions  of  compressible 
viscous  flows.  Journal  of  Comp.  Phys.,  Vol.lOl,  pp.104-129 


Mathematical  modelling  of  diagenetic  processes  in 
sedimentary  basins 

Astrid  Holstad 

Januaxy  10,  1995 


la  the  field  of  geology  diagenesis  refers  to  the  sum  total  of  processes  that  bring  about 
changes  in  a  sediment  or  sedimentary  rock  subsequent  to  deposition  in  water.  These  pro¬ 
cesses  may  be  physical,  chemical  and/or  biological,  see  [1\ 

In  this  paper  we  present  the  matnsmatical  model  and  numerical  methods  used  to  simu¬ 
late  a  subset  of  the  chemical  diagcntic  processes  within  a  sedimentary  basin.  The  processes 
are  the  dissolution  and  precipitation  or  minerals  and  decay  of  organic  carbon  present  in 
the  porous  sediments  in  contact  fTith  a  moving  fluid  phase  in  the  pores.  The  mineral  re¬ 
actions  and  the  organic  carbon  decay  processes  are  time  dependent.  In  addition  there  are 
instantaneous  equilibrium  reactions  in  the  nuid  phase. 

The  results  of  the  mineral  dissointion/precipitation  reactions  and  the  decay  of  organic 
carbon  are  changes  in  the  porosity  and  the  permeability  in  the  sediments  and  hence  the 
fluid  flow  properties  within  a  sedimentary  basin.  The  main  results  we  therefore  want  from 
our  model  are  information  about  the  behaviour  of  the  mineral  and  organic  carbon  reactions 
combined  with  information  about  how  these  reactions  affects  the  fluid  flow  properties  within 
a  sedimentary'  basin  through  geological  time. 

In  our  model  we  assume  the  sedimentary  basin  is  two-dimensional  and  represent  it 
in  space  by  qnadrilateral  elements.  Each  element  is  assigned  a  temperature,  a  mineraJ 
composition  and  a  fluid  phase  composition.  Between  the  elements  there  are  fluid  flov/  and 
dirrusion  of  ions  and  complexes  dissolved  in  the  fluid  phase. 

The  diffusion  depends  on  temperature  and  concentration  gradients.  Electrostatic  poten¬ 
tial  IS  not  considered.  Paxameters  that  influence  on  the  fluid  flow  are  viscosity,  permeabilitv 
and  pressure.  The  chemical  reactions  taking  place  in  the  basin  are  temperature  dependent. 
We  assume  Neumann  conditions  for  the  unknowns  at  the  basin  boundaries.  Typical  simu¬ 
lation  periods  are  ]  -  10  million  years.  The  temperature  and-the  fluid  flow  history  may  be 
the  result  from  a  basin  simulator. 

in  our  model  we  are  studying  the  time  dependent  development  of  the  concentration  of 
12  minerals  in  the  solid  phase  and  18  spedes  in  the  fluid  phase  per  basin  element.  The 
minerals  are  modelled  by  ordinary  differential  equations.  The  species  in  the  fluid  phase 
are  modelled  by  10  mass  balances  which  are  algebraic  equations  and  8  mass  action  laws. 
To  reduce  the  number  of  equations  to  integrate,  the  8  mass  actions  laws  are  eliminated. 
A  tymical  .simulation  i.s  10  -  100  elements,  i.e.  220  -  2200  equation.s  to  integrate  and  300  - 
3000  concentrations  to  compute. 
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The  fdl  matheiratical  model  is  thns  a  Differential  Algebraic  Equation  Svstem  fDAP 
system);  “  ‘ 

y'  =  (1,) 

“  =  s(i,r.2),  ,ii,, 

With  t  G  IR,  y  G  z  G  f  :  B  X  IR’^  g  ;  IR  x  B"  ->  and 

n  =  TiM-  rip,.  Here  n  is  the  total  number  of  unknowns,  um  is  the  number  of  mineral 
unknowns  and  is  the  number  of  fluid  phase  unknowns,  y  is  the  mineral  unknowns  and 
2  is  the  fluid  phase  unknowns. 

This  system  of  equations  has  the  following  properties; 

•  The  equations  are  very  non-linear. 

•  The  unknowns  are  tightly  coupled  and  can  inry  many  orders  of  magnitude  during  a 
Simulation  period  both  incividually  and  compared  with  each  other. 

•  mineral  reactions  switch  between  the  state  oi  dissolution,  equilibrium  and  pre¬ 
cipitation.  This  switching  is  extremely  sensitive  for  changes  in  xhe  fluid  composition. 

•  The  equations  are  very  stim 

•  It  is  not  knOoH  if  there  exists  a  unique  solution  for  the  set  of  parameters  v/e  are 
considering. 


To  solve  this  complicated  DAE-system  we  have  constructed  a  solver  based  on  the 
RADAUB  code,  see  [2].  Our  main  modifications  in  this  code  are  non-negativity  handling 
and  introduction  of  sparse  Jacobian  and  iterative  linear  solver.  Because  of  the  complexity 
of  the  equations  we  have  to  use  an  analytical  Jacobian. 

"We  will  present  numerical  results  from  the  model,  in  both  single  element  and  multi¬ 
element  basin  configurations  to  show  the  validity  of  the  model.  The  input  data  to  the  model 
are  real  data  from  sedimentary  basins  in  the  North  Sea.  The  results  will  be  presented  as 
the  time-dependent  development  of  the  mineral  composition  of  the  sediments  and  the  uore 
fluid  compo.sitiou. 
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We  introduce  the  medial  axis  transform  as  a  tool  in  the  analysis  of  geometric  structure 
of  void  space  in  porous  media.  The  medial  axis  of  any  volume  filling  object  is  intuitively 
the  “spine”  of  the  object  running  along  its  geometric  middle.  It  defines  the  inherent  shape 
of  the  object,  for  example,  the  medial  axis  of  a  cylinder  is  the  axis  of  rotational  symmetry, 
for  a  sphere  it  is  the  center  point.  The  medial  axis  for  any  n-dimensional  object  can  be 
found  by  an  algorithm  equi\-alent  to  the  following  process.  Imagine  burning  the  object  by 
igniting  a  fire  simultaneously  everywhere  on  its  entire  surface.  Let  this  fire  burn  into  the 
object  at  a  constant  rate,  reducing  it  layer-by-layer.  The  set  of  all  points  within  the  object 
where  the  fire  extinguishes  itself  is  the  media  axis.  The  fire  is  said  to  extinguish  itself  at  a 
point  if  two  components  of  the  fire,  traveling  in  opposite  directions,  simultaneously  arrive 
at  the  point. 

Volume  information  can  also  be  associated  with  each  elementary  segment  of  the  medial 
axis.  As  the  rate  of  burn  is  constant,  the  time  at  which  the  fire  burns  any  point  in  the 
object  also  gives  the  radius  of  the  largest  L2-norm  ball,  centered  on  the  point,  which  just 
fits  completely  inside  the  object.  Thus  each  point  in  the  object  is  assigned  a  unique  radius 
value,  which  we  refer  to  as  the  burn  radius;  in  particular,  each  medial  axis  point  is  assigned 
a  burn  radius,  which  we  refer  to  as  the  medial  axis  radius  at  the  given  point. 

The  medial  axis  of  an  object  gives  basic  information  concerning  its  topology  and 
geometry.  The  attraction  of  the  medial  axis  is  that  it  has  lower  dimensionality  than  the 
object  itself;  for  a  two-dimensional  object,  the  medial  axis  will  consist  of  a  union  of  lines 
and  points,  for  a  three-dimensional  object:  surfaces,  lines  and  points. 

We  apply  this  analysis  to  the  void  spaces  in  high  resolution  (5  micron),  three  dimen¬ 
sional,  tomographic  images  obtained  using  the  National  Synchrotron  Light  Source  facility 
at  Brookhaven  National  Laboratory.  Three  samples  are  analyzed,  consisting  of  two  drill 
core  samples,  Berea  sandstone  and  Danish  chalk,  and  a  sample  of  packed,  100  micron 
diameter  glass  beads.  A  Berea  sandstone  sample  was  chosen  as  it  is  a  reference  stan¬ 
dard  rock  for  petrology.  The  Danish  chalk  core  sample  is  representative  of  a  low  porosity 
rock  that  is  difficult  to  work  with  using  invasive  laboratory  techniques,  such  as  microsand- 
ing/microslicing.  The  packed  glass  bead  sample  provides  a  realization  of  a  non-overlapping 
sphere  model,  a  model  which  has  been  extensively  studied  both  theoretically  and  experi¬ 
mentally. 

In  this  poster,  we  describe  a  discrete  version  of  the  medial  axis  “burn”  algorithm 
applicable  to  digitized  tomographic  images  consisting  of  a  cubic  array  of  voxels.  Example 
medial  axes  obtained  for  the  void  spaces  in  the  data  samples  are  shown. 
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We  present  four  analyses  of  the  medial  axes  constructed  for  these  three  data  sets. 
The  first  analysis  we  present  is  for  the  distribution,  n(r),  of  burn  radii,  which  is  a  discrete 
analog  of  the  pore-size  distribution.  We  show  that  for  the  two  drill  core  samples,  the  burn 
radius  distribution  is  very  well  fit  by  exponential,  n(r)  ~  exp(r/ro).  This  conclusion  is 
also  true  for  the  distribution  of  the  subset  of  medial  axis  radii.  This  contrasts  with  the 
burn  radius  distribution  of  the  packed  glass  bead  sample  which  is  very  well  fit  by  a  normal 
distribution,  n(r)  ~  exp((r/ro)^). 

The  second  analysis  concerns  the  size  distribution  of  disconnected  medial  axis  seg¬ 
ments,  which  correspond  to  disconnected  void  volumes.  All  samples  indicate  a  composite 
distribution  consisting  of  a  power-law  distribution  of  the  smallest  medial  axis  fragments, 
followed  by  an  essentially  constant  distribution  of  medial  axis  fragments  of  intermediate 
sizes.  Finally,  both  the  Berea  and  glass  bead  data  sets  contain  a  single,  very  large  inter¬ 
connected  medial  axis  segment,  occupying  approximately  40%  of  the  pore  volume.  Such  a 
large  segment  is  absent  in  the  lower  porosity  chalk  sample. 

The  third  analysis  concerns  the  tortuosity  distribution  of  the  medial  axis  between 
parallel  planes  in  the  sample.  Assuming  that  the  medial  axis  between  two  parallel  planes 
consists  of  N  voxels,  there  are  on  the  order  of  N\  connecting  paths  to  be  considered.  For 
the  Berea  and  Bead  samples  considered  here,  the  value  of  N  can  exceed  100,000.  To 
simplify  this  CPU  intensive  computation,  only  shortest  path  connections  along  the  medial 
axis  between  the  planes  were  considered.  For  the  Berea  and  chalk  samples,  the  tortuosity 
distribution  appears  to  be  best  fit  by  an  incomplete  Gamma  distribution,  which  has  two 
free  parameters.  The  distribution  for  the  Bead  sample  does  not  appear  to  be  fit  well  by 
the  same  model. 

Finally,  considering  the  medial  axis  radius  as  a  continuous  function  along  the  medial 
axis,  we  classify  segments  of  the  medial  axis  as  throats  (local  radius  minima)  and  nodal 
pores  (local  radius  maxima).  We  present  throat  and  nodal  pore  size  distributions  and  the 
connectivity  between  them. 
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Etude  numerique  de  la  methode  de  recuperation  secondaire  du  petrole 


Mazen  SAAD^ 

On  s’interesse  a  letude  numerique  de  la  methode  de  recuperation  secondaire  du 
petrole  dans  un  gisement  constitue  d'un  seul  type  de  roche  caracterise  par  la  porosite, 
la  permeabihte  intrinseque,  les  pressions  capillaires  et  les  permeabilites  relatives. 
Cette  methode  consiste  a  injecter  un  fluide  mouillant  I’eau  dans  les  puits  d'injection 
pour  deplacer  le  petrole  vers  les  pmts  de  production. 

On  considere  que  la  phase  eau  (w)  et  la  phase  huile  (o)  sont  incompressibles 
tandis  que  la  phase  gaz  (g)  est  compressible.  La  modelisation  de  ce  type  d'ecoulement 
est  decrite  a  I'echelle  macroscopique  par 

(1)  (pdtS^+VT.'^Vw(s)+div{G^(s,Pg)g}+  v^^divVT+div(D^(s)'Vs)  =  0 

(2)  (pdiSo+VTyv^  +  diu{G^{s,Pg)g)+  v^divV^ +diviD^is)V  s)  =0 

{3)(p  dt  (Pg Sg) + Vr  .V  ipg  Vg) +div{pg  Gg  is,  Pg)g)+Pg  VgdivVT+div  ( Pg  Dg  (s) Vs ) = 0 

La  vitesse  totale  Yp  est  donnee  par 

(4)  VT  =  -KMVpg  +Gis,Pg)g+DTi8)Vs 

On  a  de  plus 

(5)  s«;+So+s^=l 

Dans  ce  systeme,  represente  la  saturation  de  la  phase  tj,  s  =  (Sj^,So),  Pg  la. 
pression  du  gaz,  v^=v^(s)  la  fraction  des  flux,  pgiPg)  la  masse  volumique  du  gaz, 
G^is,Pg)  les  termes  gravitationnels,  D.^is)  les  termes  dus  aux  pressions  capillaires 
ir=w,o,g,  T),  Mis)  la  mobilite  totale. 

Dans  le  domaine  O=]0,  L[  x  ]0,  /[,  les  conditions  aux  limites  sont 

-  sur  la  paroi  d'injection  Tg 

s^it,x)  =  l,Pgit,x)=PgOU  Vt-TI^-Q  surTg 

-  sur  la  paroi  de  production  Fg 

Pgit,x)=p,  <Pg-,iD^is)V  s).n  =  0  sur  Fg 

-  sur  la  paroi  impermeable  Fj- 

Vj  .71  =0 ;  ( D j,(s)  V  s  ).n  =  0  sur 

Le  systeme  (l)-(3)  est  un  ststeme  parabolique  degenere  ;  la  degenerescence  est  due 
au  fait  que  la  fraction  de  flux  s'annule  avec  .  Le  systeme  (l)-(3)  est  forme  de  deux 
parties  modelisant  deux  phenomenes  physiques  differents.  Le  premier  correspond  au 
transport  sous  Taction  de  la  seule  force  exterieure  d’injection  de  Teau  a  une  vitesse  ou 
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une  pression  donnee)  et  sous  I'action  des  termes  gravitationnels.  Le  deuxieme 
correspond  a  la  diffusion  sous  I'effet  des  termes  capillaires.  Afin  de  prendre  en  compte, 
aussi  precisement  que  possible,  ces  deux  phenomenes  une  methode  de  pas 
ff  actionnaires  est  utilisee. 

Dans  la  formulation  du  systeme  (l)-(3),  I'operateur  est  scinde  en  deux  parties  ;  ainsi 
la  methode  des  pas  fractionnaires  employee  ici  est  une  methode  a  deux  pas.  Nous 
decrivons  schematiquement  cette  methode. 

Au  premier  pas,  nous  donnons  une  caracterisation  pour  que  la  premiere  partie 
(correspondant  au  transport)  soit  hyperbolique  ;  ceci  donne  une  fafon  naturelle  de 
classer  les  modeles  et  permet  de  mettre  en  place  des  schemas  quasi  d’ordre  deux. 

Plus  precisement,  on  se  limitera  aux  modeles  de  permeabilites  relatives  donnant 
lieu  a  un  systeme  hyperbolique  (I'etude  de  I'hyperbolicite  de  ces  ecoulements  a  ete  es- 
sentiellement  faite  dans  [3],  [5]).  Pour  resoudre  ce  probleme,  nous  mettons  en  place 
un  solveur  approche  de  type  Roe.  Pour  ce,  nous  explicitons  une  construction  de  la 
matrice  de  Roe  pour  une  certaine  classe  de  permeabilites  relatives. 

Si  les  fliix  sont  des  polynomes  de  degre  inferieur  ou  egal  a  deux,  Vila  [6]  construit 
explicitement  une  matrice  de  Roe.  Nous  avons  generalise  ce  resultat  lorsque  les  flux 
s  ecrivent  comme  quotient  de  deiix  polynomes. 

Le  schema  de  Roe  est  d'ordre  un  et  a  trois  points.  L'inconvenient  d'un  tel  schema  est 
sa  diffusion  numerique  ;  de  nombreuses  methodes  peuvent  etre  utilisees  pour  essayer 
de  la  reduire.  Par  example,  les  limiteurs  de  pentes,  les  limiteurs  de  flux,  ou  des 
schemas  nmneriques  d'ordre  superieur  (Harten  [2]).  La  construction  de  cette  nouvelle 
famille  de  matrice  de  Roe  a  permis,  entre  autres,  la  mise  en  place  du  schema  de 
Harten.  Nous  presentons  tout  d'abord,  I'avantage  numerique  d'un  tel  schema  pour 
ecoulement  monodimensionnel,  ensuite  son  extension  aux  ecoulements 
bidimensionnels . 

Au  deuxieme  pas,  nous  traitons  les  operateurs  d'ordre  deux  en  espace  par  un 
schema  implicite  en  pression  et  semi-explicite  en  saturations.  Dans  cette  etape,  on 
etablit  une  equation  discrete  en  pression,  non  degeneree,  qui  preserve  la  consistance 
du  schema  avec  I'operateur  induit. 

Les  resultats  obtenus  sont  satisfaisants,  on  retrouve  le  comportement  attendu  : 
etalement  des  fronts  du  a  la  diffusion  capillaire,  freinage  des  fluides  du  a 
I'augmentation  de  la  mobilite  massique  du  gaz  en  compressible.  On  verifle  egalement 
la  conservation  de  la  masse  de  chaque  phase. 
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One  of  the  major  applications  of  modelling  of  flow  in  porous  media  can  be  found  in  the  petroleum  reservoir 
simulation.  However,  the  reservoir  heterogeneities,  such  as  absolute  permeability,  generated  on  fine  grids  by 
geological/geostatistical  models,  can  not  be  directly  used  for  flow  simulation,  due  to  computational  cost  and  memory 
storage  capacity.  Therefore,  efficient  techniques  are  needed  to  scale  the  fine  grid  petrophysical  parameters  up  to  the  coarse 
grid*'^. 


It  is  well  known  that  the  equivalent  gridblock  permeability  obtained  by  the  scaling-up  procedure  from  fine 
gridblocks  depends  on  the  boundary  conditions’’^.  In  the  literature,  most  boundary  conditions  considered  represent  a 
"linear"  flow  pattern.  Using  directly  these  scaling-up  results  in  reservoir  simulation  might  be  erroneous,  because  the  near 
well  behavior  can  not  be  approximated  by  the  "linear"  flow. 

On  a  reservoir  field  scale,  the  flow  region  can  be  divided  into  two  types:  a  "radial"  flow  region  with  a  high 
pressure  gradient  and  a  "linear"  flow  region  with  a  low  pressure  gradient"’.  The  "radial"  flow  region  is  usually  more 
important  for  the  prediction  of  production  forecast,  because  it  is  directly  related  to  the  well.  A  specific  scaling-up  procedure 
is  needed  in  the  near  well  region  by  imposing  a  "radial"  flow  condition.  The  transmissibility  will  be  scaled-up  in  stead  of 
the  permeability  because  it  is  more  accurate  for  flow  simulation^.  In  the  well  vicinity,  the  scaled-up  parameters  should  not 
only  consist  of  the  equivalent  transmissibility  (permeability),  but  also  numerical  productivity  index  (PI).  More  precisely, 
the  equivalent  transmissibility  (permeability)  and  numerical  PI  is  calculated  as  follow"’: 

I.  According  to  Darcy's  law,  the  equivalent  transmissibility  T^q  (Keq)  between  two  coarse  grids  is  defined  by  the 
quotient  of  equivalent  flow  Fgq,  which  is  the  sum  of  flows  at  the  fine  grid  interfaces  composing  the  coarse  grid 
interface,  and  the  difference  in  coarse  grid  pressures  Pco,  which  is  defined  as  the  storage-weighted  mean  of  the 
pressures  of  all  the  fine  grids  composing  the  coarse  grid: 

T^eq  -  ^eq^  ^Pco 

2.  By  the  definition  of  numerical  productivity  index,  the  PI  on  a  coarse  grid  is  calculated  by  the  following  formula: 

P^co  ~  P^fine  (Pfine  '  Pw)  ^  ( P^'^  "  Pw)  ~  (Pco  '  Pw) 

where  PI  fine  is  the  numerical  PI  on  fine  grid,  Pw  is  the  wellbore  pressure,  Pco  and  Pfine  are  the  coarse  grid  and  fine 
grid  wellblock  pressures,  and  Q  is  the  well  flow  rate. 

The  goal  of  the  scaling-up  procedure  is  to  reconstruct  the  fine  grid  simulation  by  a  coarse  grid  simulation, 
especially  for  the  well  production  forecast.  We  give  an  example  to  show  the  need  of  the  specific  scaling-up  procedure  in 
the  near  well  region. 

Example: 

A  heterogeneous  configuration  is  generated  on  99x99  fine  gridblocks  with  a  short  correlation  length  (Figure  la). 
Fourteen  wells  are  distributed  in  this  field  with  five  injection  wells  (11-15)  and  nine  production  wells  (P1-P9)  (Figure  lb). 
Pressure  conditions  are  imposed  at  all  the  wellbores  to  investigate  flow  rate  calculations.  A  fine  grid  simulation  is  done 
with  a  standard  finite-difference  reservoir  simulator,  where  the  numerical  productivity  indices  are  calculated  by  a 
conventional  formula^.  The  results  of  this  simulation  are  considered  as4;he  reference  solution.  A  33x33  coarse  grid  is 
constructed  based  on  the  fine  grid,  and  each  coarse  grid  includes  nine  fine  ones.  Two  scaling-up  procedures  are  used  and  the 
results  of  flow  simulation  on  coarse  grids  are  compared  to  the  reference  solution. 

The  first  scaling-up  procedure  is  the  standard  one.  i.e.,  no  flow  boundary  conditions  on  the  edges  and  constant 
pressure  at  the  inlet  and  outlet  faces^.  The  second  one  is  new:  it  includes  the  standard  procedure  for  "linear"  flow  pattern 
and  the  specific  procedure  in  the  well  vicinity  for  "radial"  flow  pattern. 

Single-phase  incompressible  flow  is  first  tested.  Figure  2  presents  the  flow  rate  for  each  well  given  by  three 
simulations:  a  fine  grid  simulation  and  two  coarse  grid  simulations  with  the  scaling-up  procedures  presented  above.  It  can 
be  seen  that  the  errors  caused  by  the  standard  scaling-up  procedure  attain  about  30%  en  the  average,  while  those  caused  by 
the  new  procedure,  which  includes  the  scaling-up  for  "radial"  flow  pattern,  are  only  about  1-2%.  This  test  shows  the  need 
to  use  the  new  procedure  when  scaling-up  in  the  vicinity  of  wells.  With  the  new  procedure,  the  error  for  well  performance 
can  be  significantly  reduced. 

Two-phase  flow  simulation  is  also  investigated.  Water  is  injected  in  an  oil  reservoir.  Figure  3  shows  the  water 
cut  for  three  production  wells  PI,  P2  and  P3.  Again,  the  simulation  with  new  scaling-up  procedure  approaches  the  fine 
grid  results  much  better.  The  breakthrough  time  can  be  calculated  accurately  by  using  the  new  scaling-up  procedure. 


It  is  clear  that  using  a  scaling-up  procedure  without  considering  the  impact  of  the  "radial"  flow  pattern  in  the  near 
well  regions  might  greatly  bias  the  flow  simulation  results  on  a  reservoir  scale.  The  simulation  accuracy  can  be 
significantly  improved  by  using  the  new  scaling-up  procedure.  This  procedure  is  very  easy  to  use  in  reservoir  simulation. 
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la).  Permeability  field 


lb).  Well  positions 


Fig.  1— Heterogeneity  and  well  configuration 


Production  well  nun±)e: 


Fig.  2-Single-phase  flow  simulation 


Fig.  3— Two-phase  flow  simulation 
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On  the  Dam  Problem  with  Leaky  Boundary  Conditions 
and  Linear  Darcy’s  Law 


We  would  like  to  present  various  results  on  the  Dam  problem  with  leaky  boundary 
conditions.  The  problem  is  to  find  a  pair  {p,x)  ^  x  L°°(Q)  such  that 

'  P  >  0,  0  <  X  <  1,  P(1  -  x)  =  0  a.e  in  Q, 

p  =  0  on  5-2, 

f  a(x)(Vp  +  xe).V^clx  -  [  f3(x,<p  -  p)4d<7(x)  <  0 

Ju  '  JSs 

e  H\a),  on  So, 

where  /?  is  a  nondecreasing  Caratheodory  function  and  a  €  yV4n(IR)). 
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Abstract 

Drying  is  an  important  unit  operation,  which  is  used  in  many  production 
processes.  An  example  of  a  branch  of  industry,  where  drying  plays  an  essen¬ 
tial  role  in  the  production  process,  is  the  clay  brick  and  tile  manufacturing 
industry. 

Shrinkage  occurs  during  drying  of  many  materials.  Moisture  concentra¬ 
tion  gradients  in  the  material  and  corresponding  gradients  in  the  amount 
of  shrinkage  will  lead  to  drying  stresses.  Controlling  these  stresses  is  im¬ 
portant  since  they  can  lead  to  deformation  or  cracking  of  the  product.  A 
correct  description  of  the  evolution  of  moisture  concentration  profiles  in  the 
material  is  complicated  by  the  influence  of  shrinkage  on  mass  transfer. 
Drying  of  clay  can  be  modelled  by  equations  expressing  mass  conservation 
of  moisture  and  solid,  together  with  a  flux  equation  for  moisture.  This  flux 
equation  is  essentially  Darcy’s  law.  where  pressure  and  permeability  depend 
on  moisture  content,  i.e.  the  quotient  of  moisture  and  solid  concentration. 
These  equations  are  supplemented  with  equations  expressing  the  deforma¬ 
tion  of  clay  due  to  gradients  in  moisture  content.  It  is  assumed  that  clay 
behaves  as  an  elastic  material,  i.e.  Hooke’s  law  applies,  where  Young’s  mod¬ 
ulus  and  Poisson’s  ratio  depend  on  the  moisture  content. 

A  careful  analysis  leads  to  a  nonlinear  parabolic  equation  for  the  moisture 
content  and  a  nonlinear  eUiptic  equation  for  the  displacements,  that  depends 
on  the  moisture  content. 

These  equations  are  discretised  in  space  by  finite  elements  and  in  time  by 
Euler’s  method  with  frozen  coefficients.  The  resulting  systems  of  sparse  lin¬ 
ear  equations  are  solved  by  preconditioned  conjugate  gradients. 

Results  of  three-dimensional  simulations  are  presented. 


The  Plane  Potential  Flow  Through  The  Inhomogeneous  Porous 

Line 

Ivan  Keglevic* 

November  25,  1994 


We  consider  the  plane  stationary  potential  flow,  i.e.  a  pair  of  functions  {v,p),v  :  — >■  IR-,  p  :  — >•  IR, 

which  satisfy 

divv(a:,y)  =  0,  (1) 

rotr(a:,y)  =  0,  (2) 

(^(a;,^/)  •  V)v(a;,T/)  +  Ivp(a;,!/)  =  0,  (3) 

where  p  is  a  positive  real  constant  (the  density  of  an  incompressible  fluid).  We  assume  that  the  flow  is  ideal 
(i.e.  the  formulae  (1)  to  (3)  are  valid)  in  the  whole  except  on  the  straight  line  y  =  0.  There,  we  have  the 
continuity  condition 

«^y+(2:>0)  =  z;j,_(a;,0),  (4) 

where  (x,  0)  is  in  porous  region  and  i'y+  and  are  the  limits  of  the  Vy  component  of  the  velocity  field  form 
the  +  (y  >  0)  and  -  (y  <  0)  part  of  the  plane.  Let  the  porosity  law  be  given  by 

7i(a;)t;j,(a;,0)  =  72(a;)(p+(a;,  0)  -  p_(2:,  0)) ,  (5) 

where  71 ,  72  :  ffi  — *•  M  are  tw’o  given  nonnegative  functions  which  do  not  vanish  simultaneously  and  and 
p_  are  the  limits  of  pressure  field  from  the  +  resp.  —  part  of  the  plane.  We  also  assume 

•i;+(a;,  0)  =  v_(a;,  0) ,  r  e  IS\supp  71  .  (6) 


If  the  velocity  in  infinity 

lim  v(x,  y)  =  Voo  =  (i'oox,  t’Dcy)  (”1) 

I(x,y)|-*oo 

and  the  functions 

71,72:1 — »[0,oo)  (8) 

are  given,  we  look  for  velocity  and  pressure  fields  v  and  p. 

The  form  (5)  of  the  porosity  law  is  derived  by  Sanchez-Palencia  [6],  Murat  [5]  and  Damlamian[1] 
and  already  used  in  WoLFERSDORF  [7]  and  in  Mikelic-Suhadolc-Veselic  [4]. 

If  the  function  72  is  zero  at  some  points,  the  hindrance  is  completely  impervious.  If  71  is  zero,  the 
hindrance  is  completely  open.  We  assume  that  71  and  72  are  piecewise  continuous  and  have  the  limit  in  the 
infinity. 

For  Voox  0  the  solution  of  this  problem  can  be  explicitely  given  (through  the  formula)  if  there  is  the 
limit 

lim  2lM  =  d<oo,  (9) 

r— >±00  72(2;) 
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<24 


(10) 


in  case  supp  ji  —M.  and  if 


lim  EM 

®->±oo  72(0;) 


=  0. 


hold  in  case  supp  71  E.  The  second  possibility  also  includes  the  cases  when  only  one  or  more  than  one 
bounded  line  segments  are  porous  or  imprevious.  The  uniqueness  of  this  solution  in  the  set  of  all  functions 
V  with  their  boundary  value  on  F  in  can  also  be  proved. 

In  the  case  Vqox  =  0  we  do  not  have  any  solutions  in  the  space  of  all  functions  v  with  boundary  values 
on  r  in  If  supp  71  consist  of  finitely  many  bounded  line  segments,  we  do  get  some  solutions  (with 

boundary  values  not  in  Ffg^(E))  but  there  is  always  more  than  one  solution,  and  all  solutions  show  the 
unexpected  property  that  the  flow  do  not  get  through  hindrance  then  round  it.  In  this  cases  the  model  do 
not  describes  properly  the  expected  physical  situation. 
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Abstract 

At  many  sites  with  subsurface  contamination,  potential  health  risks  make  remediation 
necessary.  As  a  rule,  remediation  of  subsurface  contamination  is  expensive  and  the  consequences 
of  failure  are  considerable.  Therefore,  a  thorough  assessment  of  remediation  candidates  in  the 
context  of  the  particular  site  is  warranted  before  a  selection  is  made.  At  the  heart  of  many  such 
assessments  is  the  prediction  of  fluid  flow  and  contaminant  transport  through  media  with 
multiple  scales  of  heterogeneity.  In  the  past  decade,  this  need  has  prompted  a  reevaluation  of  the 
existing  theory  for  flow  and  transport  in  porous  media.  The  central  difficulty  in  devising  new 
models  from  existing  theory  is  that  model  parameters  often  depend  on  the  scales  at  which  their 
measurements  are  made.  For  example,  an  apparent  scale  dependency  arises  as  a  result  of  the 
macroscopic  Fickian  description  for  dispersion,  which  is  evidently  a  non-Fickian  process.  Many 
non-Fickian  dispersion  models  currently  proposed  are  non-local,  wherein  the  rate  of  change  of 
concentration  in  an  elemental  volume  depends  not  only  on  a  concentration  gradient  in  the 
element,  as  for  a  local  continuum  law,  but  on  the  history  of  the  gradient  in  a  finite  region  around 
the  element.  Accordingly,  the  non-local  form  of  the  equation  is  represented  in  terms  of 
convolution  integrals.  The  non-local  form  contains  a  length  scale  that  permits  the  introduction  of 
scale  dependence. 

A  theoretical  approach  is  presented  based  on  representing  a  discrete  porous  medium  in  terms 
of  continuous  variables.  A  non-local  governing  equation,  similar  in  form  to  those  presented  by 
others  [1,  2, 3]  results.  Regardless  of  the  origin  of  the  non-local  equation,  the  convolution  integral 
may  be  approximated  discretely  as  a  network.  This  network  is  distinguishable  from  traditional 
network  models  in  three  ways;  (1)  the  flow  paths  represent  the  larger-scale  primary  advection 
paths  in  the  porous  material,  not  individual  pores,  (2)  variations  in  advection  paths  capture 
dispersion  whereas  spatial  averaging  reproduces  Darcy-scale  properties,  and  (3)  because  the 
network  is  a  discrete  form  of  the  non-local  model,  each  location  in  the  domain  may  be  connected 
to  all  other  locations,  resulting  in  unusually  high  coordination  numbers  and  a  very  dense  pattern 
of  flow  paths.  This  approximation  is  valid  at  any  scale  from  the  pore  scale  to  the  field  scale.  A 
non-local  advection  equation  will  necessarily  display  dispersion  not  found  in  its  local  continuum 
counterpart.  The  observed  scale  effect  in  apparent  dispersivity  is  reproduced  by  including 
sufficient  spatial  range  and  structure  in  the  kernel  function  of  the  convolution  integral. 
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A  discrete  network  model  was  constructed  to  numerically  simulate  non-local  advection.  Flow 
paths  are  laid  out  such  that  spatial  averages  honor  the  Darcy-scale  flow  parameters.  Variations 
in  small-scale  flow  direction  and  resistance  endow  the  network  with  the  desired  dispersion 
characteristics.  Flow  paths  are  discretized  along  their  length  to  capture  small-scale  details  and 
reduce  numerical  dispersion.  Several  notable  results  are  observed:  (1)  the  dispersion 
characteristics  of  statistically  homogeneous  media  depend  on  the  range  of  the  kernel  function,  (2) 
by  incorporating  random  variability  in  the  throat  parameters,  early  scaling  behavior  is  matched, 
(3)  the  dilemma  of  coupling  is  made  clear  by  the  network  model;  in  the  physical  system, 
processes  are  coupled  at  the  finest  scale,  a  feature  lost  by  continuum  models  but  only 
approximated  by  coarse  network  models,  (4)  model  parameters  in  these  equations  do  not  depend 
on  scale. 

The  network  represents  an  alternative  discrete  system  to  the  real  porous  medium.  Different 
network  systems  yield  the  same  Darcy-level  averages  and,  therefore,  produce  spatially  averaged 
behavior  indistinguishable  from  one  another.  In  this  sense,  the  network  is  simply  a  coarsened 
discrete  system  that  yields  the  same  dispersive  behavior  as  the  real  medium,  provided  both  are 
viewed  through  the  same  broad  band  filter. 

Research  to  this  point  shows  that,  in  principle,  a  network  may  be  constructed  that  gives  the 
best  representation  of  the  real  medium/or  a  given  level  of  resolution.  The  key  problem  for  practical 
applications  is  determining  procedures  to  calibrate  these  networks,  especially  when  non-linear 
processes  are  involved.  Efforts  are  underway  to  calibrate  the  network  model  to  specific 
laboratory  and  field  data  and  to  extend  the  approach  to  other  processes  such  as  nonconservative 
solute  transport  and  immiscible  flow. 
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Finding  effective  parameters  to  describe  a  flow  in  large  scale  heterogeneous  porous  media 
is  a  problem  of  major  concern  to  the  fields  of  hydrogeology,  petroleum  engineering  and  the 
chemical  process  industry.  In  such  systems,  the  velocity  field  u  of  the  flow  is  related  to  the 
gradient  of  pressure  Vp  by  Darcy’s  law 

u  =  —K(r)  Vp  (1) 

where  K(r)  is  the  spatially  varying  permeability  which  varies  on  length  scales  much  longer 
than  the  characteristic  pore  scale.  In  practice,  it  is  desired  to  describe  the  heterogeneous 
medium  in  terms  of  an  effective  homogeneous  medium,  i.e.: 

<  u  >= -Ke/f  <  Vp  >  (2) 

where  is  the  constant  effective  permeability  tensor. 

The  problem  of  computing  the  effective  permeability  of  a  heterogenous  medium  has 
been  the  subject  of  numerous  analytical  and  numerical  efforts  [1].  Most  of  these  works  have 
focused  on  systems  for  which  the  logarithm  of  the  permeability  has  a  normal  probability 
distribution  function.  In  our  work  [2],  we  derive  the  effective  permeability  tensor  Kg/f  of 
strongly  heterogeneous  media  for  arbitrary  distributions  of  the  log-permeability  M[K(r)].  We 
perform  a  perturbation  expansion  of  Darcy's  law  in  the  variance  <7^  of  /n[/v(r)]  and  show  how 
to  solve  the  problem  at  higher  orders  using  spectral  methods.  The  only  assumption  is  that 
the  spatially  varying  permeability  /c(r)  is  a  stationary  random  function  of  position.  This 
approach  expresses  the  effective  permeability  in  terms  of  the  moments  of  the  distribution  of 
ln[K{r)],  i.e.  Kg//  can  formally  be  computed  for  any  given  distribution  of  the  fluctuations 


of  the  log-permeability.  We  then  analyze  the  computed  effective  permeability  tensor  for 
different  particular  situations.  As  a  special  case  of  the  theory,  we  examine  Kg//  for  a 
normal  distribution  function  of  /n[K(r)]  for  both  isotropic  and  anisotropic  media.  In  this 
case,  all  odd  moments  of  the  distribution  function  are  zero.  In  the  case  of  isotropic  systems, 
a  conjecture  has  been  made  in  the  past  [3]  according  to  which,  for  three-dimensional  porous 
media,  the  scalar  effective  permeability  /Cg//  =  Kq  exp[(j^/6]  where  Kq  is  the  geometric 
mean  of  the  log-permeability.  It  is  shown  here  that  this  conjecture  is  incorrect  as  the  (re¬ 
order  term  of  Kg//  contains  additional  terms  than  those  corresponding  to  the  development  of 
the  above  formula.  Moreover,  these  additional  terms  depend  on  the  shape  of  the  two-point 
correlation  function  of  /nj/v].  The  resulting  Kg//  lies  below  the  exponential  formula.  For 
anisotropic  systems,  O(cr^)  corrections  to  the  effective  permeability  tensor  are  given  for  an 
arbitrary  orientation  of  stratification. 
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Background  and  Motivation 

More  and  more  geological  models  are  available  to  describe  the  internal  structure  of  oil  and  gas  reservoirs. 
These  models  are  the  results  of  geoscientific  work  aiming  to  integrate  the  data  and  knowledge  about  the  field. 
Generally,  these  models  are  represented  on  a  very  high  resolution  grid  (HR  Grid).  It  is  not  unusual  to  obtain  a  grid 
with  millions  of  cells.  However,  it  is  not  possible  to  solve  over  this  HR  Grid  while  fluid  flow  simulating,  because  of 
the  hugeness  of  the  linear  system  to  be  solved.  So,  one  must  upscale  petrophysical  parameters. 

The  conventional  method  aims  to  obtain  a  lower  resolution  grid  (LR  Grid)  by  coarsening  the  mesh.  The 
motivation  of  this  coarsening  is  generally  to  perform  fluid  flow  simulations  at  a  reasonable  cost.  The  results  are 
averaged  phase  pressures  and  saturations,  and,  for  compositional  modelling,  components  of  the  oil  or  gas.  These 
classical  methods  prevent  the  dynamical  upscaling  to  the  extent  that  one  homogeneizes  the  petrophysical 
parameters  before  doing  the  fluid  flow  simulation.  Some  attempts  have  been  made  to  overcome  this  difficulty  [1]. 

To  overcome  this  drawback,  an  original  method  based  on  two  grids  is  proposed.  The  principle  of  this  one 
consists  in  solving  the  parabolical  pressure  equation  over  a  LR  Grid,  and  the  hyperbolical  transport  equations 
over  a  HR  Grid.  So,  a  specific  discretisation  in  space  for  each  unknown  allows  to  keep  the  information  on  the 
distribution  of  the  saturation  while  computing  accurate  averaged  parameters  for  the  pressure  equation.  This 
method  can  be  seen  as  a  way  to  dynamically  update  the  pseudo-functions  for  the  relative  and  capillary  pressure 
curves. 


The  Model  Problem 

Let  us  consider  the  following  dead-oil  problem  with  the  hypothesis  of  incompressibility  of  the  fluids  and 
the  rock.  The  diffusion  and  capillary  pressure  are  neglected. 

The  problem  is:  find  the  pressure  p  and  the  saturation  s  such  that: 


(1) 


-div  km(s)gradp  j  -  -i- 


as 

(2)  (j)— -div 

at 


=  k 

k-^2radp 
u.,. 


) 


Pressure  Equation 


Saturation  Equation 


with  initial  conditions  in  saturation,  and  with  boundary  conditions  in  pressure  and  flow  rate,  such 
that  the  problem  is  well  posed. 

In  the  proposed  method,  each  unknown  will  have  its  own  discretisation.  An  IMPES  scheme  is  used.  So, 
we  apply  an  implicit  scheme  for  the  pressure  equationon  on  the  LR  Grid,  and  an  explicit  scheme  for  the  saturation 
equation  on  the  HR  Grid. 


The  Dual  Mesh  Algorithm 

Let  us  describe  an  extension  of  the  Dual  Mesh  Algorithm 
presented  for  the  pressure  equation  [2]. 

We  consider  the  update  from  to 

H  is  relative  to  the  LR  Grid  and  h  to  the  HR  Grid  (cf  fig  1 ) 

Step  1 :  Homogeneization  of  the  petrophysical 
parameters  km(s[])  over  the  LR  Grid 

Step  2:  Resolution  of  the  pressure  equation; 
calculation  of  pj)*' 


Reservoir 

"  ^  HR  Grid 
LR  Grid 

Dual  Mesh  Representation 


-div  km(s"  )grad(pH''' )  =  fw  +  f, 


fig  1 


3^ 


step  3: 


Reconstruction  of  the  pressure  and  the  flow  rate  thanks  to  and 


Step  4  :  Calculation  of  the  saturation  by  resolution  of  the  transport  equation  over  the  HR  Grid: 


=  sfl  +  (j)At"div 


gradph""’ 


+  L 


Step  3  -  Reconstruction  of  the  flow  rate  and  the  pressure 

Let  us  describe  into  details  the  step  3.  This  is  one  of  the 
most  important  in  this  algorithm. 

In  order  to  take  into  account  the  heterogeneities  of  the 
porous  media,  we  propose  an  original  method  by  solving  for  each 
cell  (called  Mj  -  cf  fig  2)  of  the  LR  Grid  a  local  problem  using  the 
same  interpolator  as  the  one  used  to  model  the  global  system. 

There  are  several  options  to  define 

-  the  region  r(Mj)  around  each  cell  Mj. 

-  the  boundary  conditions  on  r(M|) 

We  choose  first  r(Mj)  =Mj 

The  problem  of  the  boundary  conditions  is  probably  the 
most  difficult  point.  After  taking  boundary  conditions  in  pressure, 
such  that  the  flow  rate  Q  is  kept  on  the  boundary  of  Mj,  we  decide 
to  take  Neuman  conditions  on  Mj  with  one  Dirichiet  condition. 


Numerical  Results 

This  Dual  Mesh  Method  was  successfully  implemented  on  synthetic  cases.  In  our  poster,  one  could  notice 
how  the  interpolator  follows  the  variation  of  the  pressure  due  to  the  discontinuity  permeability  field.  At  last,  the  Dual 
Mesh  Method  is  completely  validated  in  homogeneous  case,  with  a  good  estimation  of  the  saturation. 


Conclusion 

The  paper  describes  an  algorithm  which  allows  to  solve  the  transport  equation  over  the  HR  Grid  after 
solving  the  pressure  equation  on  a  LR  Grid,  and  obtaining  the  flow  rate  at  a  smaller  scale  by  solving  local  problems. 
This  algorithm  is  implemented  for  a  dead-oil  problem.  These  results  are  encouraging  and  could  be  extended  to  more 
general  multiphase  flow. 


Nomenclature 


fw,  fo: 

Symbols  pertaining  to  the  LR  grid  are  written 

source  terms  vector 

in  capital  letters,  those  for  the  HR  grid  are  written  in  small  letters. 

p\  interpolated  pressure  vector  (over  the  HR  Grid) 

k-. 

permeability  tensor 

p: 

pressure  vector  (over  the  HR  grid) 

krw: 

relative  permeability  of  the  water 

s: 

water  saturation 

mj: 

current  cell 

P-  w; 

viscosity 

m: 

mobility 

4)  : 

porosity 
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Contribution  a  I'etude  du  comportement  de  la  solution  d'equations  non  lineaires 

de  diffusion-convection 

J.B.Betbeder,  L.Ldvi,  A.Plouvier,  G.Vallet 

Laboratoire  de  Math6inatiques  Appliqu6es,  URA  1204  CNRS.  IPRA,  avenue  de  I'Universitd,  64000  PAU. 


La  moddlisation  mathdmatique  des  phdnomSnes  naturels  rencontrds  en  mdcanique  des  milieux 
continus  fait  appel  k  la  traduction  de  lois  de  conservation  (de  masse,  de  I'dnergie,  ...).  On  est  alors 
amend  k  considdrer  des  dquations  non  lindaires  de  type  divergentiel,  dites  dquations  de  continuitd. 
EUes  sont  gdndriquement  formuldes,  au  sens  des  distributions  dans  Q  =  ]0,T[xQ  domaine 

bornd  rdgulier  de  , 

I)  soit  sous  la  forme  d'un  probldme  de  Cauchy  de  diffusion-convection 

(I)  ^^-A(j)(U)  +  Div  4^(U)G  =f, 
at  L  - 

P(U)(0,.)  =  P(Uo_), 

auquel  on  adjoint  des  conditions  aux  limites  approprides, 

(p ,  (f)  et  'F  sont  au  moins  des  fonctions  continues,  (j)  dtant  croissante  et  p  strictement  croissante) 

II)  soit  sous  la  forme  d'une  dquation  hyperbolique  non  lindaire  du  premier  ordre  ; 

3U 

(H)  — -f-Div  T(U)G  =f, 

at 

U(0,.)  =  Uo, 

lorsque  I'effet  de  transport  est  preponddrant  par  rapport  k  I'effet  de  diffusion-dispersion.  Ceci  est  par 
exemple  le  cas  dans  une  phase  fluide  hdtdrogdne,  lorsque  Ton  ddsire  suivre  I'dvolution  de  la  fraction 
massique  d'un  constituant.  En  effet,  la  diffusion-dispersion  moldculaire  de  ce  dernier  (due  k 
I'existence  de  gradients  de  concentration  au  sein  de  la  phase)  peut  etre  "ndghgeable"  au  regard  de  la 
Vitesse  de  filtration  du  fluide. 

(la  fonction  T est  gdndralement  supposde  de  classe  C^,  Uq  bomde  et  k  variation  bomde) 

Lorsque  I'on  cherche  ^  montrer  I'unicitd  de  la  solution  de  tels  probldmes  par  des  mdthodes  de 
troncature  dans  V,  il  est  ndcessaire  de  prendre  en  compte  les  6ventuelles  ondes  de  chocs  d'une 
solution. 


1.  i)  Ainsi,  lorsque  la  fonction  (j)  est  injective  (probl&me  de  type  parabolique  degenere),  s'il 

existe  une  solution  U  de  I'dquation  (1)  telle  que  — ^  €  L'  (£2) ,  alors  cette  solution  est  umque.  Bien 

ot 

sur,  on  n'est  jamais  assure  de  I'existence  d'une  telle  solution  "forte".  Cependant,  moyennant  certaines 
hypotheses  de  rdgularitd  sur  la  donnee  initiale  et  sur  les  non-lindarites  p  et  (|),  G.Vallet  a  su  dormer 
dans  [12]  un  rdsultat  d'existence  d'une  solution  "forte"  pour  I'dquation  intdgro-diffdrentielle 

ap(u) 

AA/TT^  _  T/'TTM": _ f  ^T/'TT/',. 


-A6(U)  =  J(U)Exp  -Jj(U(s,.))ds  dansQ, 


vdrifiant  des  conditions  de  bord  de  type  Fourier.  Sa  demarche  s'appuie  sur  im  resultat  ddveloppe  dans 
[2]  etmontre  mSme  que  cette  solution  est  enfait  dldment  de  W^’”  (0,T;L'  (H)). 


1.  ii)  Lorsque  1'  on  n'est  pas  en  mesure  d'dtabhr  I'existence  d'une  solution  "forte"  de  (I),  il  est 
ndanmoins  possible  de  montrer  I'unicitd  d'une  solution  U  vdrifiant  :  (j)(U)  e  H’ (Q)  et 
P(U)  e  BV(0,T;L' (Q))  nL“  (Q).  La  ddmonstration  repose  sur  la  validation  de  la  formule  de 
ddrivation  de  la  composition  d'applications  dans  I'ensemble  BV  des  fonctions  k  variation  bomde.  Pour 


cela,  on  montre  que  I'onde  de  chocs  (au  sens  de  la  limite  approximative)  est  n6gligeable  par  rapport  h 

la  mesure  de  Radon  (cf.  par  exemple  [6],  [3]  ou  encore  [10]). 

3t 

Enfin,  I'adaptation  des  resultats  expos6s  dans  [4],  a  permis  a  G.Gagneux  &  M.Madaune-Tort 
[7]  de  montrer  qu'rme  solution  faible  de  I'dquation  (I)  pos6e  dans  Hj,  (^2)  v6rifie  implicitement  une 
formulation  de  Kruzkov.  Sans  hypotheses  supplementaires  sur  les  non-lmearit6s,  son  unicitd  est  alors 
prouv6e  en  reprenant  I'idee  de  la  methode  gen6ralement  ddvellopde  dans  le  cas  d'xme  Equation 
hyperbolique  non  lineaires  du  premier  ordre.  (voir  aussi  [5]  pour  des  probl&mes  de  barrages,  [11]  lors 
que  le  tenseur  de  diffusivite  depend  de  revolution  du  probl^me  ou  encore  [12]  dans  le  cadre 
d'^quations  de  la  thermodynamique) 

II.  La  solution  de  I'equation  (11)  est  obtenue  par  I'dtude  du  comportement,  lorsque  £  tend  vers 
0"^,  de  la  suite  des  solutions  des  problemes  (P^ )  de  "diffusion  lente" : 

9Up  r  ^ 

(Ig )  -  eA(t)(UE )  +  Div  )  G  =  f , 

ot  L  . 

Ue(0,.)  =  Uo. 

associes  k  des  conditions  de  bord  Dirichlet  sur  tout  ou  une  partie  de  la  frontibre. 

Reprenant  I'id^e  de  la  methode  de  viscosity  artificielle,  on  est  amend  k  estimer  la  solution  du 

probleme  (P^)  dans  I'espace  (Q)  nL”  (Q),  dans  les  cas  de  solutions  "fortes",  voire  dans 

I'espace  BV(Q)  n  L°°  (Q),  dans  le  cas  de  solutions  faibles.  Ainsi,  en  raisonnant  par  compacitd,  on 
dtablit  I'exitence  d'une  solution  faible  entropique  U  a.  variation  bomee.  (cette  technique  de 
perturbations  singulidres,  introduite  par  [1],  a  notamment  dtd  reprise  par  M.J.Jasor  [8]  pour  des 
probldmes  k  effet  de  puits  et  par  L.Levi  dans  le  cas  de  probldmes  unhatdraux  [9]). 

RemarQue  :  Puisque  I'onde  de  chocs  (au  sens  de  la  limite  approximative)  n'est  pas  ndgliaeable  par 

3u  au 

rapport  aux  mesures  de  Radon  —  et - ,  i=l..n,  on  comprend  mieux  la  ndcessitd  d'une  formulation 

at  ax; 

entropique,  qui  permet  de  "controler"  les  sauts  de  U  au  niveau  d'lme  onde  de  chocs  (relation  de 
Rankine-Hugoniot)  et  d'assurer  ainsi  son  unicitd. 


References 


[1]  C.Bardos,  A.Y.LeRoux.  J.C.Nedelec  :  Pint  Order  Quasilinear  Equations  with  Boundary  Conditions,  Comm,  in  Partial 
Differential  Equations,  4  (9),  1017-1034.  (1979). 

[2]  Ph.Benilan,  R.Gariepy  ;  Strong  solutions  in  LI  of  Degenerate  ParaboUc  Equations,  J.  Diff.  Equa.,  a  paraitre. 

[3]  J.B.Betbeder,  G.Vallet  :  La  question  de  I'unicite  pour  une  classe  de  problemes  devolution  non  lineaires  integro- 
diff&entiels,  C.R.Acad.Sci.Paris,  T.317,  Serie  1.  1993,  pp,319-322. 

[4]  J.Carrillo  :  1.  Unicite  des  solutions  de  type  Kruskov  pour  des  problemes  elliptiques  avec  des  termes  de  transport  non 
lineaires,  C,R,Acad.  Sci.  Paris,  T.303,  Serie  1,  n'5,  1986.  pp.  189-192. 

[5] J.CaiTillo  :  On  the  Uniqueness  of  the  Solution  of  the  Evolution  Dam  Problem,  Nonlinear  Analysis.  Vol.  22.  n^5,  1994, 
pp,573-607. 

[6]  G.Gagneux,  M.Madaune-Tort ;  Sur  la  question  de  I'unicite  pour  des  inequations  des  milieux  poreux,  C.R.Acad.Sci.  Paris, 
T,314,  Serie  1,  1992,  pp.605-608. 

[7]  G.Gagneux,  M.Madaune-Tort  :  Unicite  des  solutions  faibles  d'equations  de  diffusion-convection,  C.R.Acad.Sci.Paris, 
T.318,  Serie  1, 1994,  pp.919-924. 

[8]  M.J.Jasor  :  Perturbations  singuKeres  d’equations  non  lineaires  de  diffusion-convection  modelisant  des  ecoulements 
diphasiques  incompressibles  en  milieu  poreux.  These  de  doctorat  d'universite,  Universite  de  Pau  (1992). 

[9]  L.Levi ;  Equations  quasi  lineaires  du  premier  ordre  avec  contrainte  unilaterale,  C.R.Acad.Sci.Paris,  T.317,  Serie  1,  1993, 
pp  1133-1136 

[10]  L.Levi  :  Perturbations  singulieres  d'une  inequation  parabolique  degeneree,  publi.  interne  93/15  URA  1204-CNRS, 
Univ.  de  Pau. 

[lljA.Plouvier  :  Sur  une  classe  de  problemes  d'evolution  quasi  lineaires  degeneres,  Revista  Matematica  de  la  universidad 
Complutense  de  Madrid,  a  paraitre. 

[12]G.'Vallet ;  Strong  Solution  of  a  Nonlinear  Integro-differential  Parabolic  Equation,  Ann.  Fac.  Sc.  Toulouse,  a  paraitre. 


Sur  la  correction  non-lineaire  de  la  loi  de  Darcy  pour 
les  ecoulements  en  milieux  poreux  a  faible  nombre  de 

Reynolds 

Mouaouia  Firdaouss,  Jean-Luc  Guermond 

1.  Introduction.  -  On  demontre  une  conjecture  annoncee  dans  Wodie-Levy  [5]  et  Mei- 
Auriault  [2]:  dans  les  milieux  poreux  periodiques  dont  les  inclusions  sont  de  taille  equivalente 
a  la  periode,  la  correction  non  lineaire  a  la  loi  de  Darcy  est  quadratique  en  fonction  du  nombre 
de  Reynolds  =  pV PqR j ou  VPq  est  Fechelle  de  gradient  de  pression  ma-.  nscopique,  I 
Fechelle  spatlale  microscopique,  p  la  masse  volumique  du  fluide  et  p.  sa  viscosite  dynamique. 
On  montre  que  I’hypothese  d’isotropie  faite  dans  [5]  et  [2]  n’est  pas  necessaire,  mais  peut  etre 
j-gjij-jplacee  par  une  liypothese  d'invanance  par  rapport  au  sens  du  flux.  Des  tests  numeriques 
illustrent  ce  result  at. 

2.  Position  du  PROBLEME.  -  On  suppose  le  milieu  poreux  periodique.  On  se  donne 

pour  cellule  periodique  un  parallepipede  de  1R“  engendre  par  les  vecteurs  (ti, . .  ■,tn)  qu’on 
suppose  lineairement  independants,  c  est  a  dire  dete(ii, . . .  ,tn)  =  mes(n)  0,  ou  e  est  la 
base  canonique  de  IR”  (ie.  en  se  donnant  e,  on  fixe  1  echelle  microscopique).  La  base 
(ti,...,tn)  n’etant  pas  forcement  orthonormee,  on  definit  la  base  covariante. 

Pour  i  =  1, . . .  n  on  note  Sfl'  la  face  de  0.  passant  par  Porigine  de  IR"  (qu’on  munit  de  sa 
structure  affine)  et  telle  que  SD'.P'  =  0;  on  pose  aussi  dCif  =  dQ,~  +  t,-.  Soit  (Iq  un  ouvert 
connexe  de  fi,  de  frontiere  reguliere  (disons  Lipscliitz  et  d  un  seul  cote  de  sa  frontiere). 
Pour  i  =  1 , . . . ,  n  on  pose  =  dO-f  (T  dflo  et  on  suppose  que  flo  est  compatible  avec  la 
periodicite,  c’est  a  dire  dflfp  =  30. t  L-  Ou  definit  enfin,  dQo^s  —  dClQ\  Uj_;i  30, -p.  On 
suppose  que  mes(30o,5)  7^  0  et  mes(30-p)  #  0  pour  i  =  1, . .  .  ,n.  On  introduit  maintenant 
le  cadre  fonctionnel  habituel:  soient  Ho  j,(no)  =  {u  G  H^(Oo), 

V  =  {v  e  Ho,p(0o),  divu  =  0}. 

Dans  un  premier  temps  on  considere  le  probleme  suivant;  pour  A  6  S'n(0,l)  (la  sphere 
unite  de  E”,  ie.  on  fixe  le  gradient  de  pression  macroscopique  de  reference),  trouver  uo(A)  G 

V  tel  que 

(2.1)  Vu  G  V,  (Vuo(A),  Vu)  =  (A,u). 

On  introduit  maintenant  le  tenseur  de  permeabilite  Kq  G  £(E”,  IR")  tel  que  /GA  =  J^^  Uo{X). 
Ce  tenseur  est  symetrique  defini  positif.  Dans  le  cadre  de  la  theorie  de  1  homogeneisation 
on  interprete  la  relation  KqX  =  uo(A)  comme  la  loi  de  Darcy  (c/.  Bensoussan-Lions- 
Papanicolaou  [1],  .Sanchez-Palencia  [3]  et  Tartar  [4]).  L  objectif  de  ce  travail  est  d  etudier  la 
correction  non  lineaire  de  cette  relation  dans  le  cadre  des  equations  de  Navier— Stokes. 

On  se  place  maintenant  en  dimension  u  —  2  ou  3.  On  introduit  la  forme  trilineaire 
b{u,v,w)  =  {u.Vv,wy,  b  est  continue  sur  Hop(Do)^,  pour  n  =  2  ou  3,  on  note  la  norme 
de  b  sur  pi^oY-  On  s’interesse  au  probleme  de  Navier-Stokes:  pour  e  >  0  et  A  G  5„(0, 1) 
trouver  iq(A)  G  V  tel  que 

(2.2)  Vu  G  P,  (Vu,(A),  Vu)  +  efi(tq(A),u,(A),u)  =  (A,u). 

Si  e  <  co(no)'V2c!,mes(flo),  ou  co(Do)  est  la  constante  de  Poincare  du  domaine  Do,  A  existe 
une  solution  unique  au  probleme  ci-dessus.  Par  la  suite,  on  caracterise  la  dependence  de 
ttj(A)  par  rapport  a  A  pour  e  petit. 
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3.  Correction  non  lineaire  de  la  loi  de  Darcy  -  On  introduit  maintenant 
rhypothese  ceatrale  de  ce  trava.il: 

(H)  On  suppose  que  pour  tout  A  G  5„(0,  1)  et  tout  e  suffisamment  petit,  ^^(A)  satisfait; 

Remarque  1  :  Cette  hypothese  signifie  qu’un  renversement  du  gradient  de  pression  macro- 
scopique  implique  un  renversement  de  la  vitesse  de  filtration  sans  alteration  notable  de 
son  module,  et  ce  independamment  du  nombre  de  Reynolds  (suppose  petit),  le  resultat  de 
I’experience  pouvant  depeudre  de  la  direction  (non  orientee)  du  gradient  de  pression  en  ques¬ 
tion. 

Remarque  2  :Si  Oq  est  invariant  par  symetrie  centrale,  (H)  est  verifiee  et  la  correction  0[e^) 
est  exactement  nulle  independamment  de  t. 

Remarque  3  :  Si  fio  fst  invariant  par  rapport  aux  symetries  Si,  i  —  l,...,n,  telles  que 
Si[ti)  =  ti  et  Si{tj)  =  —tj  si  j  ^  i,  alors  Qq  est  invariant  par  symetrie  centrale. 

Le  resultat  central  de  ce  travail  est  le  suivant 
Theoreme  3.1  Si  rhypothese  (H)  est  satisfaite,  on  a 

(3.1)  VA  G  -5,1(0, 1),  f  iq(A)  =  AqA  -f  0{e^). 

J  r^o 

Ce  resultat  signifie  que  pour  e  petit,  la  vitesse  de  filtration  verifie  quasiment  la  loi  de  Darcy 
avec  le  tenseur  de  permeabilite  Kq,  la  norme  de  I’erreur  etant  bornee  par  ce^,  ou  la  constante 
c  ne  depend  que  de  la  geometrie  du  domaine. 

On  pent  s’afFranchir  de  I’hypothese  (H)  si  recoulement  de  filtration  est  monodimension- 
nel.  Plus  precisement,  on  se  place  dans  le  cas  ker(Ao)'^  =  lm(/{^o)  c’est  a  dire  Kq 

se  reduit  a  la  multiplication  par  une  constante;  on  montre  alors 

Theoreme  3.2  Dans  le  cas  monodimensionnel ,  on  a 

(3.2)  VA  G  {A/|C|,-^i/|iii),  [  u,(A)  =  AToA  +  Cl(e^). 

4.  Illustrations  NUMERIQUES  -  Les  resultats  annonces  ici  font  en  fait  suite  a  une  serie 
d’experimentations  numeriques  en  dimemsion  2  sur  des  reseaux  periodiques.  Nous  rapportons 
ci-dessous  quelques  unes  de  ces  experiences  afin  d’illustrer  les  tlieoremes  3.1  et  3.2. 

Les  figures  1,  2,  3  et  4  resument  des  resultats  numeriques  significatifs  pour  quatres  con¬ 
figurations  differentes.  A  gauche  on  represente  revolution  de  |  u^(A)  —  A’cA]  normalise  par 

son  maximum  en  fonction  du  nombre  de  Reynolds  relatif  Re/ Rmax  pour  differentes  valeurs 
du  Reynolds  maximum  Rmax-  A  droite  des  figures  on  a  represente  la  configuration  periodique 
correspondante.  Dans  tous  les  cas  on  a  choisi  |ii|  =  1. 

La  figure  1  concerne  un  reseau  rectangulaire  de  rapport  |ii|/|t2|  =  2.  Le  tenseur  Kq  n’est 
pas  spherique,  le  milieu  homogeneise  n’est  done  pas  “isotrope”  pour  les  equations  de  Stokes 
(il  n’entre  pas  dans  le  cadre,  de  [2],  [5]).  Ce  reseau  est  a  symetrie  centrale;  il  verifie  done 
(H).  On  verifie  sur  la  figure  que  pour  R,  <  20,  la  correction  non  lineaire  est  parfaitement 
quadratique.  Le  cas  presente  correspond  a  A  =  e'^^A^ 
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La  figure  2  coiicerne  un  reseau  rectangulaire  qui  n’est  pas  a  symetrie  ceiitrale  inais  qui 
satisfait  (H)  experinaentalement.  On  verifie  encore  id  que  la  correction  non  lineaire  est 
c^uadratique  (on  a  pris  A  = 

La  figure  3  concerne  un  cas  monodimensionnel.  II  n’y  a  pas  de  reversibilite  (re.  (H)  n’est 
pas  verifiee),  mais  le  theoreme  3.2  s’applique.  Independamment  du  sens  de  I’ecoulement  et 
pour  /?£  <  10  on  verifie  que  la  correction  non  lineaire  est  quadratique. 

La  figure  4  concerne  un  cas  bidimensionnel  fabrique  de  telle  sorte  qu’il  ne  satisfait  pas 
(H).  Ce  cas  a  ete  assez  difficile  a  exhiber.  Pour  A  ==  on  verifie  que  la  correction  non 

lineaire  est  lineaire  en  fonction  du  nombre  de  Reynolds. 
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Figure  1:  Cas  anisotrope  a  symetrie  centrale 
Figure  1:  Anisotropic  case  with  central  symmetry 


Figure  2:  Cas  anisotrope  avec  (H)  verifiee  experimentalement 
Figure  2:  Anisotropic  case  for  which  (H)  is  satisfied  experimentally 
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Figure  3:  Cas  monociimensionnel  sans  reversibilite  du  flux 
Figure  3:  Monodimensional  case  without  reversibility  of  the  flux 


Figure  4:  Cas  anisotrope  avec  (H)  non  verifiee 
Figure  4:  Anisotropic  case  for  which  (H)  is  not  satisfied 
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Abstract 


This  paper  considers  the  Cauchy  problem  for  the  following  nonlinear  system 


Vt  -  Ux  =  0 

(1.1)  Ut  -  (J(v,s)x  +  au=  0 

T 

with  bounded  measurable  initial  data 

(1.2)  (v,u,s)/t=o  =  (vo(x),uo(x),so(x))  , 

where  a,  (3,  x  are  nonnegative  constants.  When  p  =  0,  system  (1.1)  can  be  used  to  model  the 
adaiabatic  gas  flow  through  porous  media,  where  v  is  specific  volume,  u  denotes  velocity,  s 
stands  for  entropy  and  a  denotes  pressure.  Its  form  in  Euler  coordinates  is  also  a  model  of 
isothermal  unsteady  two  phase  flow  in  pipelines.  In  this  paper  we  study  the  global  generalized 

solution  for  the  case  P  =  0. 

When  P  0,  which,  when  written  in  Euler  coordinates,  can  be  used  to  model  the 
chemically  reacting  flow,  where  again  v  is  specific  volume,  u  denotes  velocity  but  s  is  the 
mass  fraction  of  one  mode  of  the  two-mode  gas  and  f(v)  is  a  given  equilibrium  distribution  in 
V.  In  this  case,  t  denotes  a  reaction  time.  We  show  that  the  solution  of  the  equilibrium  system 

(1.3)  vt-ux  =  0 

ut  -  a(v,f(v))x  -t-  au=  0 

is  given  by  the  limit  of  the  solutions  of  the  viscous  approximation 
vt  -  Ux  =  evxx 

( 1 .4)  Ut  -  a(v,s)x  +  au=  e  Uxx 


as  the  dissipation  and  the  reaction  time  x  go  to  zero. 


MODELLING  OF  A  HORIZONTAL  AND  A  VERTICAL- 
FRACTURED  WELL 


B.Cvetkovic,  G-  Halvorsen,  E.  L0w,  Rele  -  Reserv'oir  Tedmology  Section 
IFE  -  Kjeller,  Norway 


A  horizontal-well  and  vertical-fractured- well  production  rate-time  analysis 
needs  a  different  approach  and  treatment.  Reservoir  simulation  data  and  a 
well-production  data  may  be  compared  to  the  modelling  data.  The 
difference  to  the  wellbore  rate-time  data  will  be  analysed  and  discussed. 
The  rate-time  and  pressure-time  comparison  data  analysis  should  provide 
more  realistic  future  rate  prediction.  Comparison  analysis  will  be  based  on 
developed  type-curves  data,  reservoir  simulation  rate-time  data  and  w'ell 
production  data.  The  rate-time  data  will  be  compared  and  obtained  results 
discussed.  The  approach  will  comprise  modelling,  simulation  analysis  of 
both  a  horizontal-well  and  a  vertically  fractured  well  in  a  single  porosity 
reservoir  that  may  be  extend  to  the  dual  porosity  case.  The  comparison 
analysis  will  start  with  a  simplified  horizontal-well  case  and  in  a  future  the 
proposed  simplified  model  can  be  extended  to  a  model  that  will  include 
heterogeneities  (vertical  to  horizontal  permeability  differences,  faults  and 
multilayers  or  composite  structures). 
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Resolution  numerique  d'un  probleme 
de  restauration  biologique  en  milieu  poreux 


Ch.H.  Bruneau^  -  P.  Fabrie^’^ 


On  s'interesse  4  la  resolution  numerique  d'un  problfeme  modelisant  la  restaura¬ 
tion  biologique  par  des  bacteries  B  d'un  milieu  poreux  pollue  par  im  substrat  miscible 
S  en  presence  d'oxygene  0  [4].  Le  modele  etudie  est  celui  de  Borden  et  Bedient  [1]  ou 
I'adsorption  des  substrats  en  solution  est  supposee  lineaire  et  instantanee,  et  les 
termes  de  reactions  qui  regissent  le  developpement  des  bacteries  et  la  degradation  du 
polluant  sont  issus  de  la  cindtique  de  Monod  [1]  : 


(1)  V((A(S)+D(V)).VS  )+y.VS=-OBs  {0,S,B) 

dt 

(2)  V((A(0)-hB(V)).V0)+y.V0=  -<I>Bo(0,S,B) 
3 1 

(3)  §  =  Rb{0,S,B) 

(4)  ^i(S)V=-K.(Vp-l5(S).e,) 

(5)  v.y=o 


La  vitesse  de  I'ecoulement  est  notee  V,  la  porosite  du  milieu  d) ,  la  permeabilite  du 
milieu  K,  R  est  le  facteur  de  retard  du  substrat  S,  p  et  P  la  viscosite  et  la  densite  du 
fluide,  A  le  tenseur  de  diffusion  moleculaire  et  D  le  tenseur  de  dispersion 


IXV)  =  ai 


uv ' 
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.  -uv 
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ou  et  sont  les  coefficients  de  dispersion  longitudinale  et  transversale. 

La  vitesse  de  filtration  V  etant  donnee  par  la  resolution  de  I’dquation  de  Darcy  (4) 
et  de  I'equation  de  continuite  (5),  les  equations  (1),  (2),  et  (3)  sont  decouplees  par  une 
methode  de  pas  fractionnaires  [3]  :  le  systeme  d'equations  differentielles  issu  des 
termes  de  reaction  est  alors  resolu  par  une  methode  de  Runge  Kutta  explicite  d'ordre 
quatre.  La  principale  difficulte  reside  dans  le  traitement  des  equations  d’advection  dif¬ 
fusion.  En  effet  lorsque  le  transport  domine  la  diffusion,  les  methodes  usuelles  ge- 
nerent  soit  de  la  diffusion  numerique  soit  de  la  dispersion  numerique  source 
d'  etalements  des  fronts  et  d'instabilites. 


1  LAMB,  Universite  Bordeaux  I  -  33405  Talence 

2  LEPT-ENSAM  URA  873 


A  partir  du  schema  de  Tackacs  [6],  qui  est  du  type  Lax-Wendroff  auquel  un  terme 
de  troisieme  ordre  en  espace  est  ajoute  minimisant  ainsi  les  erreurs  de  dispersion  et  de 
dissipation,  on  construit  par  la  technique  des  limiteurs  de  flux  de  Roe  ou  Davis  [5],  un 
schema  explicite  quasi  d'ordre  deux  a  variation  totale  d^croissante  pour  I'approxi- 
mation  des  termes  de  transport  [4].  Les  termes  de  diffusion  sont  discretises  de  fa?on 
impHcite  par  un  schema  centre.  Les  resultats  obtenus  sont  trds  satisfaisants. 

L'extension  bidimensionnelle  du  schema  de  Tackacs  est  basee  sur  la  discretisation 
par  la  methode  des  volumes  finis  sur  un  maillage  cartdsien  de  I'dquation 

(6)  iC_^y(2)*(V)VC)+V.VC=  0  ;  I>*(V)  = 

dt  2 

issue  du  developpement  de  Taylor  a  I'ordre  deux  en  temps  de  la  solution  C  de  I'equa- 
tion  de  transport 

(7)  — +  V.VC=0. 

dt 

Le  flux  D*(  V)VC  est  discretise  et  limite  de  fafon  a  avoir  un  schema  L°°  stable. 

Des  essais  numeriques  ont  ete  effectues  pour  le  probleme  modele  (7)  avec  des 
champs  de  vitesse  k  divergence  nulle  pour  lesquels  on  connait  une  solution 
analytique.  La  comparaison  de  la  solution  calculee  avec  la  solution  exacte  montre  que 
le  schema  est  faiblement  diffusif  et  minimise  les  effets  dispersifs  dus  au  maillage  tout 
en  conservant  la  masse.  C'est  ce  schema  qui  est  propose  pour  resoudre  la  partie 
transport  du  systeme  (1)  -  (5). 
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ASYMPTOTIC  BEHAVIOUR  OF  SOLUTIONS  OF  NEUMANN 
PROBLEMS  FOR  NONLINEAR  NONVARIATIONAL  ELLIPTIC 
EQUATIONS  IN  DOMAINS  WITH  COMPOSITE  STRUCTURE 


Alexander  A. Kovalevsky 

Institute  of  Applied  Mathematics  and  Mechanics 
Donetsk,  Ukraine 


The  asymptotic  behaviour  of  '"‘-generalized  solutions  of  Neumann  problems  for 
elliptic  equations 

-  Y,  ^  ai{x  u) a{x  .u)  =  f  ,  x  e  fla  ,  5  =  1,2,..., 
i=l 

is  studied.  The  following  qualitatively  different  kinds  of  domains  Vlg  C.  0.  are  consid¬ 
ered: 

1)  strongly  connected  perforated  domains; 

2)  weakly  connected  domains  of  framework  type  with  thin  channels; 

3)  perforated  domains  with  so  called  accumulators. 

The  first  kind  of  domains  0.,  is  characterized  with  existence  of  a  sequence  of  linear 
continuous  extension  operators  p,  :  — >  IF*’"'!!])  such  that  sup^  ||ps||  <  oo. 

In  this  case  solutions  of  the  problems  under  consideration  converge  to  solution  of  Neu¬ 
mann  problem  for  an  equation  of  the  form 

n 

—  ^  ^  (7,(x.  Vfz  i  -r  c\n{x.u)  =  af  ,  x  6  Q  .  ft  >  0  . 

In  the  second  case  domains  fl,  have  representation  Oj  =  mes  ffj  — >  0 

and  there  exist  sequences  of  linear  continuous  extension  operators  ■p^'p  :  IF*’"‘(fl^'’)  — + 
— such  that  supj  lips’ll  <  cc  (z  =  1,2).  In  this  situation  solutions  of  the 
initial  Neumann  problems  converge  to  solution  of  an  equation  with  operator  defined  on 
(lF''"‘(n))2. 

In  the  third  case  domains  Oj  have  representation  'fla  =  ,  mes  H,  0, 

there  exists  a  sequence  of  operators  p'P  :  TF^'"'!!!^”)  — »■  lF*''"(fl)  such  as  in  the  second 
case,  open  sets  Es  are  unions  of  nonintersected  domains  {j  G  J^),  which  are  called 
as  accmnulators.  In  this  case  solutions  of  the  problems  under  consideration  converge 
in  a  certain  sense  to  solution  (u.  c>)  G  lF^''"(fl)  x  of  a  problem  for  functional- 

differential  system  of  the  form 

— ftt  5(x,  u  —  ly )  -r  ft2[n(:r,  tl')  —  f]  =  0  .  x  G  , 

n 

—  ai{x.  Vu)  -f  ftoa(x.  u )  -h  ftoulx.  b’)  =  (fto  +  ft?)/  ,  x  G 

(ft,-  >0.  i  =  0.1.2)-. 


Stockage  de  gaz  reactif  dans  un  milieu  poreux 
Modelisation  mathematique  et  numerique 

Victor  Duval,  Hamid  Ghidouche  &  Claude  Basdevant 
URA  742,  C.N.R.S.,  Universite  Paris  Nord,  France 


Le  stockage  rapide  de  grandes  quantites  de  tritium,  isotope  radioactif  et  toxique  de 
I’hydrogene,  est  un  probleme  important  pour  les  installations  utilisees  dans  la  fusion  ther- 
monucleaire  controlee.  line  des  methodes  preconisee  par  les  ingenieurs  du  CEA  consiste  a 
faire  absorber  ce  gaz  par  une  poudre  metallique,  cette  reaction  etant  reversible.  L’absorption 
du  gaz  se  fait  par  une  reaction  d'hydruration  qui  est  exothermique.  Cette  reaction  se  pro- 
duit  sous  certaines  conditions  de  temperature  et  de  pression.  Durant  Thydruration  la  poudre 
gonfle,  modifiant  la  porosite  du  milieu.  Le  stockage  maximum,  c’est  a  dire  Thydruration  max¬ 
imale,  demande  de  controler,  par  un  s}-steme  de  refroidissement  inclus  dans  les  parois.  la 
temperature  du  four  de  stockage:  en  particulier  pour  empecher  la  formation  en  entree  de 
four  d’un  bouchon  d’hydrure  bloquant  la  progression  du  gaz.  Le  processus  de  destockage 
est  symetrique,  a  un  leger  hysteresis  pres  .  La  reaction  de  deshydruration  est  activee  par 
chauffage  de  la  poudre.  Le  phenomene  est  done  caracterise  par  des  couplages  non-lineaires 
raides  entre  dynamique  et  chimie  et  entre  chimie  et  thermique. 

La  modelisation  mathematique  et  numerique  du  probleme  a  ete  entreprise  dans  le  but 
d’optimiser  la  forme  du  four  et  le  controle  de  la  reaction.  Elle  est  basee  sur  I’equation  de  la 
cinetique  chimique  de  la  reaction  d'hydruration/deshydruration,  sur  le  bilan  de  masse  et  sur 
le  bilan  d'energie.  Si  le  probleme  est  relativement  standard  dans  le  cas  du  transport  d'un  seul 
gaz,  il  est  moins  classique  quand  on  prend  en  compte  un  melange  tritium  -  helium.  L'helium 
qui  peut  apparaitre  a  des  concentrations  non  negligeables  dans  le  systeme  est  inerte  dans 
la  poudre.  Nous  proposons  dans  ce  cas  une  modelisation  de  la  dynamique  du  melange  en 
definissant  des  vitesses  propres  pour  chacun  des  gaz  ;  ces  vitesses  sont  calculees  elles-memes 
a  partir  d'une  unique  loi  de  Darcy  fondee  sur  une  viscosite  dynamique  moyenne. 

La  modelisation  num&dque  est  faite  en  volumes  finis  avec  des  schemas  conservatifs.  semi- 
implicites  en  temps  et  decentres  amont  en  espace.  Des  resultats  d’existence  et  de  stabilite  ont 
pu  etre  etablis  ainsi  qu'un  principe  du  maximun  discret  analogue  du  principe  du  maximum 
obtenu  sur  le  modele  mathematique  continu. 

Le  code  de  calcul,  prenant  en  compte  la  geometric  detaillee  du  prototype  industriel 
(poudre,  masses  metalliques  du  four  et  volumes  d’entree  et  de  sortie)  donne  des  resultats  en 
accord  avec  les  experiences. 
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MACROSCOPIC  PERMEABILITY  OF  THE  SYSTEM  OF  THIN  FISSURES 
FILLED  BY  MATERIAL  WITH  RANDOM  PERMEABILITY  TENSOR 

PANASENKO  G.  P. 

Consider  the  simliest  model  of  lattice-structures:  two-dimensional  rectangular 
lattice. 

Definition  1.  The  union 

=  U^T-oo  (  eB?  \  \x2-ke  \  <  eii/2  } 

U  {{xi,X2)  \\xi-k£\  <  £\l/2  }  ) 

is  called  the  two-dimensional  rectangular  lattice. 

Thus  the  rectangular  lattice  is  a  union  of  thin  strips  of  the  width  s/x  stretched 
in  each  coordinate  direction  and  forming  the  e-  periodic  system  in  each  dimension. 

We  also  denote 

2:2)  e  I  I  xz-j  -ke\  <  e/x/2  } 
the  unions  of  horisontal  (j  =  1)  and  vertical  {j  =  2)  strips,  so 

Be,y.  =  Bl^^  U 

Let  C  be  a  domain  with  the  bondarj’  dG  €  which  is  independent  of  e  and 

IX. 

Definition  2.LetUe,^^{x)  is  a  sequence  of  functions  from  L2{Be,yLC[G)  ,  uo{x)  e 
L2{G).  One  sais  that  L  -  converges  to  uq{x)  at  f\G  if  and  only  if 

^ymeas{Bs,^  fl  G) 

The  normahsation  factor  1/ y^meas{Bs,^  H  G)  is  necessary  because 
iil|lL2(B£  nG)  =  ^/meas{B.,^  n  G).  Notice  that  L-convergence  is  not  a  con¬ 
vergence  in  common  sense  because  the  domain  depends  on  small  parameters. 

Consider  a  lattice  B.^^^  and  for  each  strip  Bl  =  {(xi,a:2)  €  |  1  xz-j  - 

ke  \  <  £/x/2  },  j  =  1, 2  we  associate  a  random- valued  constant  (2  x  2)— matrix  , 
independent  of  e  and  /x.  All  A^j  are  independent  in  aggregate  and  have  the  same 
discrete  distribution: 

P[A,, 

where 

S=1 


=  =  Ps,  s  =  l,...,r, 

=  L  aW  =  (A^*))^  >  0. 


Let  n  be  an  intersection  of  the  strips: 

and  be  a  fixed  constant  matrix  such  that  >  0.  We  pose 

A  =  on  the  set  II  and  A  =  A^j  in  each  strip  -B^\n  without  11. 

Let  G  be  a  domain  with  the  bondary  dG  G  C°°  which  is  independent  of  e  and 
f  E  C^{G).  Consider  the  problem  with  random  coefficients 

(1)  —div  (  A  grad  Ue,^)  =  f{x)  , 

(2)  (  A  grad  ,  n)  =0,  for  x  E  n  G, 


(3)  =  0,  for  X  E  n  dG. 

Problems  (l)-(3)  simulate  a  problem  of  permeability  of  a  fissured  rock  filled 
with  porous  substance  ,  with  A  being  the  permeability  tensor  of  the  substance  in  the 
fissures  ,  Ue^^j,  is  a  microscopic  pressure,  and  uq  is  a  macroscopic  pressure.  Numerical 
solution  of  problems  (l)-(3),  with  c<<l,  g.  «  1,  is  very  difficult  since  the  step 
size  of  the  grid  must  have  an  order  much  less  than  £.  The  realization  of  the  standard 
homogenization  procedure  is  also  impeded,  since  the  problem  on  a  cell  depends  on 
the  small  parameter  g,  and  in  order  to  solve  it  numerically,  we  must  select  the  step 
size  of  the  grid  to  be  much  less  than  g.  Hence,  an  asymptotic  investigation  of  the 
problem  is  needed. 

Introduce  the  notation  A  =  (dij), 

r 

an  =  ^P50.5(aii  —  ai2a^2^a2i),  di2  =  0, 

Sr=l 

r 

U22  =  ^Ps0.5(a22  —  a2lUii^Ul2),  U21  =  0. 

5=1 


Let  uo  be  a  solution  of  the  homogenized  averaged  problem 
(4)  —div  (  A  grad  uq)  =  f(x)  ,  x  E  G,  uolsG  =  0, 

Theorem  For  each  5  E  (0, 1)  the  estimate  takes  place 


^meas{Bs^^  fl  G) 


>  +  y/gf^  ^  C£^ . 


d  0^ 


Application  of  convective  transfers  on  geology 

B.  Gerard,  J.J.  Royer,  C.  Le  Carlier  de  Veslud 

CRPG,  15  rue  Nd  des  Pauvres.  BP  20,54501  Vandoeuvre-Les-Nancy,  France. 
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Convective  heat  transfer  in  porous  media  is  of  fundamental  importance  to  a 
number  of  geological  applications,  such  as  oil  recovery,  water  supply  management 
in  hydrogeology,  geothermal  exploitation,  ground  heat  storage,  radioactive  waste 
management,  ground  water  flow  modelling  and  is  also  of  interest  in  environmental 
sciences  and  geophysics.  But  problems  encountered  in  that  kind  of  application 
are  different  from  those  encountered  in  laboratories  for  several  reasons;  i)  large 
spatial  extension,  ii)  3D  complex  geometry,  iii)  scarce,  irregular  or  partial  data,  iv) 
geometry  of  objects  not  known.  \’)  large  amount  of  data.  In  such  situation,  there 
is  no  uniqueness  of  solution  and  consequently,  robust  benchmarks  are  necessary  to 
select  the  most  appropriated  model  compatible  with  the  available  data. 

The  purpose  of  this  work  is  to  present,  firstly,  a  global  process  for  3D  modelling 
of  transfer  problems  encountered  in  geology,  and  secondly,  indirect  procedure  based 
on  direct  measurements  for  validating  the  model.  In  a  first  step,  the  geometry  is 
modelled  using  surfaces  defining  the  limits  of  the  studied  object  (layers,  volume,...) 
and  the  discontinuities  (fractures,  faults).  Then,  properties  of  the  objects  (perme¬ 
ability.  thermal  conductivity,..)  are  interpolated  into  the  space  from  a  given  set 
of  ’’control  points”.  This  step  is  done  by  the  GOCAD  software,  developed  at  the 
Ecole  de  Gologie  de  Nancy,  a  powerful  tool  designed  for  interactive  3D  modelling  of 
surfaces  and  volumes,  of  complex  geological  structures. 

In  a  second  step,  a  special  scheme  allowing  modelling  of  coupled  heat  and  mass 
transfer  in  2D  and  axisymmetric.  anisotropic  and  heterogeneous  porous  media  has 
been  developed.  The  media  is  saturated  with  a  single  fluid  phase.  Darcy  flow  char¬ 
acteristics  are  assumed  for  the  liquid  phase.  Solid  characteristics  depend  on  the 
temperature.  A  dimensionless  formulation  is  used  to  simplify  the  heat  and  Darcy 
equations.  For  two  dimensional  and  axisymmetric  problems,  the  resulting  equations 
are  strictly  similar  and  can  be  solved  numerically  using  the  same  procedure.  Nu¬ 
merical  studies  are  in  agreement  with  results  obtained  by  theoretical  approach  for 
homogeneous  isotropic  medium. 

This  approach  has  been  applied  to  the  site  of  Soultz-sous-Forts,  selected  by 
the  European  Hot  Dry  (HDR)  geothermal  project  as  a  pilot  zone  for  exploiting 
low  enthalpy  energy.  It  is  situated  east  of  the  Merkwiller-Pechelbronn  oil  field  , 
along  the  western  limit  of  the  upper  Rhine  Graben  in  northern  Alsace  (France). 
Geophysical  and  geochemical  studies  carried  out  on  Rhine  Graben  show  evidences 
of  water  circulation,  within  the  Triassic  and  altered  granitic  rocks.  Firstly,  the 
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sedimentary  cover  is  asymmetric  with  a  depth  of  1400m  along  the  west  border  of 
the  graben  compared  to  3500m  on  the  eastern  part,  causing  a  general  circulation 
pattern  from  east  to  west,  associated  with  a  higher  surface  heat  flow  on  the  west 
(up  to  ISOmW  compared  to  70-80mW  on  the  east).  Secondly,  the  vertical 
heat  flow  decreases  vertically  from  a  mean  constant  value  of  150  mW  m~'^  in  the 
sedimentary  cover  to  a  lower  value  of  70  and  30  mW  at  1700-2000m  depth,  in 
the  bedrock  in  boreholes  GPKl  and  EPSl,  respectively.  Finally,  fluid  circulation  is 
presently  observed  in  the  granitic  fractured  bedrock.  Regional  geochemical  studies 
suggest  that  these  fluids  result  from  a  mixing  of  low  salinity  surface  water  from  the 
west  w'ith  higher  salinity  deep  formation  water  from  the  east.  The  decrease  of  the 
geothermal  heat  flow  at  depth  is  likely  attributed  to  fluid  circulation  through  the 
sediment  cover.  Physical  parameters  and  boundary  conditions  are  sometime  difficult 
to  define.  Thermal  conductivity  has  been  measured  directly  on  a  conductivimeter 
from  rock  samples  whilst  permeability  is  estimated  from  empirical  approximations 
and  injection  tests  in  boreholes.  Boundary  conditions  are  defined  from  hydrological 
and  geophysical  information. 

To  obtain  a  ’’good  model",  several  tries  have  been  necessary,  implying  the  modi¬ 
fication  of  the  geometry,  the  physical  parameters  and  the  boundary  conditions.  The 
final  model  shows  three  hydrothermal  circulation  systems.  The  first  one  is  expressed 
at  the  surface  next  to  Baden-Baden.  The  second  has  a  complex  path,  consisting  of 
six  convective  cells  extending  from  the  Rhine  axis  to  Pechelbronn.  The  third  and 
smallest  one  consists  of  two  cells,  which  originate  in  the  western  part  of  the  graben 
and  continue  to  Hochwald.  The  maximum  filtration  velocity  observed  is  about 
lm.yr“h  The  temperature  obtained  from  the  numerical  model  are  in  good  agree¬ 
ment  with  those  from  direct  borehole  observations  and  with  temperature  estimated 
from  fluid  inclusions,  suggesting  that  the  temperature  field  would  be  stable  from 
the  formation  of  the  fluid  inclusions  to  the  present  time.  The  geochemical  study 
of  water  shows  that  fluid  exchange  has  been  taken  place  betw'een  the  Triassic  and 
the  fractured  granite  layers,  resulting  from  a  mixing  of  high  salinity  formation  fluids 
and  surface  water.  These  two  last  arguments  are  in  favour  of  convective  circulation 
and  show  that  water  e.xchange  between  the  Triassic  and  the  fractured  granite  layers 
has  been  taken  place.  The  C'’’  fluid  dating  provides  ages  ranging  from  5000  to  20 
000  yr.,  compatible  with  the  velocity  calculated  by  model. 

The  above  methodology  shows  that  fluid  transfer  in  natural  porous  media  implies 
a  feedback  procedure  between  the  modelling  and  the  available  data.  The  solution  is 
not  unique,  but  indirect  methods  (fluid  dating,  etc...)  can  be  used  to  constrain  the 
final  model.  Such  procedures  are  used  to  provide  temperature  estimation  at  depth 
for  geothermal  purpose.  Together  with  geometrical  reconstitution  throughout  times, 
they  give  also  interesting  information  for  oil  formation  and  migration. 


AAo 


The  Explicit  Solution  of  a  Free  Boundary  Problem  for  a  Nonlinear 
Absorption  Model  of  Mixed  Saturated-Unsaturated  Flow 


A.  C.  BRIOZZO-D.  A.  TARZIA 

Departamento  de  Matematica,  F.C.E.,  Universidad  Austral 
Paraguay  1950,  (2000)  Rosario,  ARGENTINA. 


ABSTRACT 

In  wetted  soils,  zones  of  saturation  develop  naturally  in  the  vicinity  of  impermeable  strata, 
surface  ponds  and  subterranean  cavities.  Hidrology  must  be  concerned  with  transient  flow  trough 
coexisting  unsaturated  and  saturated  zones.  Models  of  advancing  saturated  zones  necessarily  involve  a 
nonlinear  free  boundary  problem. 

A  closed-form  analytic  solution  is  presented  for  a  nonlinear  diffusion  model  under  conditions 

of  ponding  at  the  surface.  The  soil  water  diffusivity  is  restrict  to  the  special  funcional  form 

D('6'i= _ ^ _  where  6  is  the  water  content  field  to  be  determined  and,  a  and  b  are  positive 

^  ^  (b-0)2 

constants.  The  explicit  solution  depends  of  a  parameter  C  (determined  by  the  data  of  the  problem), 
according  to  two  cases  :  1  <  C  <  or  C  >  Cj  ,  where  Cj  is  a  constant  which  is  obtained  as  the 
unique  solution  of  an  equation.  This  result  complements  the  study  given  in  P.  Broadbridge,  Solution 
of  a  Nonlinear  Absorption  Model  of  Mixed  Satured— Unsatured  Flow”,  Water  Resources  Research, 
26(1990),  2435-2443. 


On  some  new  models  of  nonhomogeneous  fluid  filtration 

in  porous  media 


V.N.  MONAKHOV 
Novosibirsk,  RUSSIA 

In  this  report  I  propose  some  new  models  of  nonhomogeneous  fluid  filtration  in 
porous  media.  The  models  based  on  Navier-Stokes  and  boundary  layer  equations. 
Existence  theorems  are  proved  as  well  as  numerical  results  are  given. 
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Multiscale  models  :  A  tool  to  describe  the  porosity 
of  cement-based  materials  and  to  predict 
their  transport  properties 
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The  porosity  of  cement-based  materials  covers  several  orders  of  magnitude,  gener¬ 
ally  ranging  from  nanometers  to  a  few  micrometers.  Therefore,  it  is  difficult  to  relate 
porosity  and  transport  properties  for  such  materials.  In  this  paper,  a  method  based 
on  the  representation  of  the  pore  structure  by  networks  of  capillary  tubes  is  proposed. 
The  basic  idea  is  to  sort  the  pore-size-distribution  into  several  classes  and  to  distribute 
each  class  on  networks  of  various  sizes.  The  different  networks  are  then  superposed 
by  following  rules  based  on  renormalization  theory.  This  model  has  been  applied  to 
correct  the  pore-size  distribution  (PSD)  obtained  from  experimental  mercury  intrusion 
porosimetry  curves.  Using  this  “corrected”  PSD,  transport  coefficients,  such  as  satu¬ 
rated  (water,  air)  permeability  and  mixed  (water  vapour-f-water)  transfer  coefficients, 
can  be  estimated  using  numerical  methods.  The  estimated  values  are  compared  with 
experimental  results. 


